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ABSTRACT 


The Pseudo Wigner-Ville Distribution is a time-frequency representation of an input 
time signal and is ideally suited for portraying non-stationary signals. A working com¬ 
puter program is presented and the effect of preprocessing and postprocessing data ma¬ 
nipulations is shown. The program has been developed for analyzing data for use in 
machinery condition monitoring and diagnostics and will be a valuable asset for ana¬ 
lyzing transient machinery. A practical example showing pump speed variations with 
time is also presented. Due to the fact that the Pseudo Wigner-Ville Distribution can 
be used to analyze both steady state and transient operations, along with the fact that 
it can be calculated on virtually any computer, this method could revolutionize machin¬ 
ery condition monitoring and diagnostics. 
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THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research may not 
have been exercised for all cases of interest. While even - effort has been made, within 
the time available, to ensure that the programs are free of computational and logic er¬ 
rors, they cannot be considered validated. Any application of these programs without 
additional verification is at the risk of the user. 
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I. INTRODUCTION 


The physical condition or state of health of machineries which operate in short du¬ 
ration cycles is not known with any degree of accuracy. Maintenance on these ma¬ 
chineries is being conducted periodically in order to avoid failures and prolong the useful 
operating life of the equipment. These machineries, since they operate for only short 
periods of time, can be characterized as transient. Additionally, machineries which are 
not operating in a steady state condition are also transient. This includes rotating ma¬ 
chinery transitioning from one speed to another. 

In order to assess the physical condition of machinery without complete disassem¬ 
bly, a physical measurement of its vibrations is conducted using an accelerometer. Other 
sensors, such as temperature or pressure transducers, could also be used. There are 
other methods, including motor current signature analysis on electrically driven ma¬ 
chinery and wear debris analysis which could be used.[Ref. 1] However, vibrations are 
used predominantly for machinery condition monitoring. The vibrations are recorded 
in the time domain. 

The time domain representations of vibrations may be decomposed into a summa¬ 
tion of sine waves and can be identified in the frequency domain (see Figure 1 on page 
2 [Ref. 2]). As long as the physical characteristics of the machinery are known, these 
sine waves or frequency components can be directly attributed to physical events oc¬ 
curring within the machinery. For example, a shaft rotating at 3600 RPM (60 Hz) with 
a 10 tooth spur gear will have a gear mesh frequency of 600 Hz (10x60 = 600). 
Therefore a frequency component of 600 Hz may be attributed to the gear. As the speed 
of the shaft changes the gear mesh frequency will also change. All physical components, 
including bearings, couplings, impellers, rotors, etc., may be related in a similar manner 
to frequencies dependent upon the shaft rotation speed. Therefore, for transient ma¬ 
chinery, as speeds vary with time the frequency components will also vary. 

There is a need for a method to represent the time dependent events which occur 
with machinery operating in transient modes. At each instant in time as the speed of the 
machinery changes the frequency content will also change. The Pseudo Wigner-Ville 
Distribution is the method which was chosen to portray these time dependent changes. 
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AMPLITUDE 


AMPLITUDE 



Figure 1. Time and Frequency Domains 

The Pseudo Wigner-Ville Distribution is a three dimensional (time, frequency, and 
amplitude) representation of an input signal and is ideally suited for portraying transient 
phenomena. The Wiener Distribution has been used in the areas of optics [Refs. 3 , 4, 
5] and speech [Refs. 6, 7], Wahl and Bolton [Ref. S] are using it to identify structure- 
borne noise components. Flandrin et. al. [Ref. 9] recently proposed its use in the area 
of machinery condition monitoring and diagnostics, while Forrester [Ref. 10] is investi¬ 
gating its use in gear fault detection. 

The Pseudo Wigner-Ville Distribution can be used to portray both transient non¬ 
stationary phenomena as well as starionary phenomena and therefore can be used for 
machinery condition monitoring of all machinery. This includes machinery operating in 
a steady state condition. Due to the time independent nature of steady state operating 
conditions, the frequency content will be constant for all time. The benefits of using the 
Pseudo Wigner-Ville Distribution and monitoring all machinery are enormous. Ma¬ 
chinery which has never before been monitored now can be. Additionally, monitoring 
now is not limited to just a given steady state operating condition, all speeds can be 
monitored, the effects of different modes of vibration can be investigated, and a more 
thorough evaluation of the machinery condition can be obtained. This all translates into 
a tremendous economic savings. Here-to-fore unmonitored machinery now can be 
monitored and machinery monitoring is not limited to just steady state operating con¬ 
ditions. 






The WIGFUN computer programs presented require that data be collected in the 
time domain and be digitized. Once it is digitized, then virtually any computer can take 
the digitized data and analyze it. The only hardware required is a transducer with a 
power supply, an analog to digital converter, and a computer. 
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II. THE PSEUDO WIGNER-VILLE DISTRIBUTION 


A. WIGNER DISTRIBUTION FUNCTION 

The Wigner Distribution Function (WDF) was First introduced by E. Wigner in 1932 
[Ref. 11]. Claasen and Mecklenbrauker, in a three part series of papers, developed 
mathematical equations for the WDF in continuous time [Ref. 12] and discrete time 
[Ref. 13]. They also showed relations with other time-frequency transformations [Ref. 
14], For the continuous time case using two different complex signals, r(r) and s(/). the 
cross-Wigner Distribution can be formed. The cross-Wigner Distribution is defined as: 

tl'DF rs (i. w) = J°° e~ jaJ y(t + f )s'(/ - f )dr (1) 

where: 

r — >■(/); a complex time signal 

s = sd): a complex time signal 

i = time 

o) - frequency 

* = complex conjugate 

The auto-Wigner distribution is defined as: 

1VDF S ,,(/, co) = f°° e~ j ^s{i + ^ )s{i - ± )dr (2) 

"-OC ^ 

Since the objective of this research is to accurately portray the time dependent frequency 
characteristics of a single input signal, only the auto-Wigner distribution will be used. 
From the frequency domain, the auto-Wigner distribution is defined as: 

V ; DF S ' S (co, i) -J~ S'Sio + j- )S*(cu - \ )d- (3) 

where: 

5 = S(u>) ; Fourier Transform of s(/) 
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There are several properties of the WDF which are important to note. As may be 
seen above, the WDF in the frequency domain and the WDF in the time domain are 
related as follows: 


lVDF SiS (co,t)=lVDF, iS (!,cD) (4) 

A time shift of a signal is a time shift of the WDF: 

lVDF^ico) = U 'DF S s (t - t, co) (5) 

A frequency shift of a signal is a frequency shift of the WDF: 

WDF e ,a, s e ,o, s {t, co) = IVDF S<s (t, co-Cl) (6) 

It follows that a time and frequency shift of a signal is both a time and frequency shift 
of the WDF: 


]VDF e^),F% {l -4 r - ") = n 'DF StS (t- T, CD - Cl) (7) 

Equation (7) is extremely important when we consider that these changes in time and 
frequency are exactly what we want to portray for our use in characterizing transient 
phenomena for machinery condition monitoring. It has been shown that the WDF can 
discriminate the frequency content of a signal at discrete times. 

Integrating the WDF over time, frequency, and both time and frequency provides 
signal energy information. The integral of the WDF over frequency at a specific time 
yields the instantaneous signal power at that time: 

J 1VDF s Jt, co)dco = I s(0 1 2 (S) 

The integral of the WDF over time at a specific frequency yields the energy density 
spectrum of a signal at that frequency: 

J"" 1VDF S ' S (t,co)dt= |S(cd)| 2 (9) 
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The integral of the WDF over the whole plane, both time and frequency, yields the total 
energy in the signal: 


°° f°° lVDF s s (t, <o)did(o = | |s| | 2 (10) 

-OO —OO 

The computer program listed in Appendix A does not accurately maintain the energy 
information in the original input signal. In the preprocessing program the input time 
signal is windowed and amplified. In the postprocessing program the distribution is 
averaged over both time and frequency. Additional properties of the WDF are listed in 
references 12 and 15. 

The discrete time auto-Wigner distribution as developed by Claasen and 
Mecklcnbrauker [Ref. 13] is defined as: 



U'DF S s (t, cu) = 2V e' J2w "s(i + t)s*(/ - r) (11) 

?= —DC 

The time and frequency shift properties of the continuous time WDF described by 
equations (5). (6), and (7) remain valid for the discrete time WDF. Yen [Ref. 15] defines 
the discrete time auto-Wigner distribution for a sampled time signal s(r), for 0< i< T, 
as: 


r -1 

1 • 4 -COT , 

IFDF,,,(/,co)=^2/ ~ *' + * ('- T ) ( 12 ) 

T=0 

Both equations (11) and (12) are similar. In either case the WDF is basically the Fourier 
transform of an auto correlation of a signal. 

B. WIGNER-VILLE DISTRIBUTION 

In 1948. Ville proposed the use of analytic signals in time- frequency representations 
of real signals [Ref. 16]. An analytic signal is a complex signal which contains both real 
and imaginary components. The advantage of using an analytic signal is that in the fre¬ 
quency domain the amplitudes of negative frequency components are zero. This satisfies 
mathematical completeness of the problem by accounting for all frequencies, yet does 
not limit our practical application since only positive frequency components have a 
practical interpretation. 
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An analytic signal may be formed from a real signal by several methods. These 
methods may be grouped into either a time domain formulation or a frequency domain 
formulation. In the time domain a typical formulation uses the Hilbert transform and 
may be expressed as: [Ref. 17] 

s{t) = s r (t)+JH{sM} (13) 


where: 

s(/) is the resulting analytic signal which is complex 
s,(t) = real component of a signal 

= imaginary component of a signal 

//{$,(')} is a Hilbert transform and is defined as: 





t* 0 


/ = 0 


In the frequency domain a typical formulation uses the Fast Fourier Transform (FFT) 
and is expressed as: [Ref. IS] 


s{i) = FfT -1 [S(a>)] 


(14) 


where: 

S*(cj) 

Siw) = 2 S R {io) 

0 

s^(w)-Ffrwo] 


cj = 0 

W- l.“- 1 

otherwise 


The distribution resulting from an analytic signal being processed through the Wigner 
distribution is commonly termed a Wigner-Ville Distribution. 


C. PSEUDO W1GNER-VILLE DISTRIBUTION 

The Wigner-Ville Distributions of most signals are very complicated and difficult to 
interpret since the input signals consist of many components. The most annoying 
characteristic of the Wigner and Wigner-Ville Distribution is the presence of cross terms. 
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Jones and Parks [Ref. 19] described the production of a cross term for the sum of two 
signals. The Wigner Distribution of the sum of two signals, r{i) + s(i), is defined as: 

WDF r+s>r+s (t, to) = tt’DF^t, «) + 2 Re[ 11 'DF,Jt, ft))] + WDF s s (t, a>) (15) 

The annoying cross term results from the cross-Wigner distribution WDF ri and is lo¬ 
cated midway between the auto terms. As the input signals consist of a summation of 
greater numbers of individual components, the number of cross terms also increases. 

There are several approaches for the removal, or deemphasis, of these cross terms. 
For machinery condition monitoring applications the presence of cross terms in the 
Wigner-Ville Distribution is not disastrous as long as they can be identified as such. 
However, they will make the resulting distributions more difficult to interpret. Usually 
an averaging process is performed in order to produce a more understandable and in¬ 
terpretable representation of the input signal. 

There are two methods for making the Wigner Distribution more presentable. 
Claasen and Mecklenbrauker describe the application of a sliding window in the time 
domain before calculating the Wigner Distribution [Ref. 12]. The resulting distribution 
is called a Pseudo Wigner Distribution. A second option is to smooth the Wigner Dis¬ 
tribution with a sliding averaging window in the time-frequency plane, sometimes re¬ 
ferred to as a Smoothed Wigner Distribution. Since the effect of this second method is 
essentially the same as the first, the removal of undesired components, the resulting 
distribution from either method will be called a Pseudo Wigner Distribution. In both 
cases the result is to deemphasize components arising from calculations and to empha¬ 
size deterministic components. Obviously, averaging a Wigner-Ville Distribution will 
result in a Pseudo Wigner-Ville Distribution. 

There has been a fair amount of research as to the optimum smoothing algorithm. 
References 12, 15, 20, 21, 22, and 23 all discuss various methods available. Due to the 
fact that the resulting Pseudo Wigner-Ville Distributions were to be used to identify time 
varying frequency components of machinery signatures, a sliding exponential window in 
the time-frequency plane is used in this work for the averaging process. This decision 
is supported by Nuttall's work which demonstrated that smoothing must be applied to 
both time and frequency [Ref. 23]. 
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III. THE W1GFUN COMPUTER PROGRAMS 


A. BRIEF DESCRIPTION 

The computer program developed in this research is included as Appendix A. The 
computer code begins with a summary of the programs, subroutines, and symbols used. 
The computer code then lists the WIGFUN Fortran programs, is followed by an 
alphnumeric listing of the subroutines, and is concluded with a listing of the data format 
conversion program referred to in Appendix B. A flow chart of the WIGFUN programs 
follows: 



Figure 2. Flow Chart of WIGFUN Computer Programs 


The program is written in Fortran 77 and is designed to be interactive so few, if any, 
modifications should be required for implementation by other users. The programs may 
not have been validated for all cases of interest. The 'include' statements in the Fortran 
programs cause the compiler to include the indicated file at compilation. The included 
files are individual subroutines for the most part, making the size of the individual pro¬ 
grams and subroutines more manageable. The plotting routines all use 
CA - D1SSPLA ® software [Ref. 24], 
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A description of the process of transferring data from a source to the Digital VAX 
computer used is provided in Appendix B. The only requirement for the data to be used 
in the computer program is that it is in the format read in subroutine 
DA TAIN. INCLUDE. 

All of the graphs presented in this thesis resulted from analog signals obtained in the 
laboratory and transferred into the VAX computer using the procedure provided in 
Appendix B. 

Figure 3 on page 11 and Figure 5 on page 13 show two Pseudo Wigner-Ville Dis¬ 
tributions. Figure 3 is the result obtained from operating on two sine waves generated 
by function generators and then added together. It can be clearly seen that the fre¬ 
quency content remains constant over all time, as expected for a sine wave. Figure 4 
on page 12 is the time signal input for Figure 3. Figure 5 is a linear chirp which was 
obtained from a sine wave whose frequency was linearly varied with a triangular wave. 
This figure demonstrates that the Pseudo Wigner-Ville Distribution is able to represent 
time varying frequencies. Figure 6 on page 14 is the time signal input for Figure 5. 
Figure 7 on page 15 shows different viewing angles of the two sine waves seen in 
Figure 3. The top left view shows all three axes, the top right view looks down the fre¬ 
quency axis, the bottom left view looks down the time axis, and the bottom right view 
looks down the amplitude axis at a single plane located at one third of the maximum 
amplitude. These four views provide a good, fast presentation of the distribution. A 
more detailed contour plot, as in Figure 8 on page 16, gives a better indication of the 
time varying characteristics of the input signal. A practical machinery monitoring ex¬ 
ample is presented in chapter 4. A description of the computer program listed in Ap¬ 
pendix A and the options available follows. 
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Figure 3. Pseudo Wigner-Ville Distribution of Two Sine Waves 























Amplitude 



Figure 4. Time Signal for Two Sine Waves 
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Figure 8. Detailed Contour Plot of Linear Chirp 
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B. LIMITING DECISIONS 


As can be seen in Figure 2 on page 9, the first Fortran program is WIGFLN1. 
WIGFL’Nl reads in the raw time data, preprocesses it, and outputs an analytic data file. 
There are several data preprocessing options available within WIGFUN1 an*, each will 
be discussed. However, prior to preprocessing, several limiting decisions must be made. 

The first decision is the number of data points to use. The maximum number of 
data points which may be processed is set by the array sizes dimensioned throughout the 
computer code. As printed the listing is set for 512 data points. If more data points are 
required then larger arrays should be dimensioned. Consequently, more computer stor¬ 
age space is required and the run time is longer. It was found that for initial applications 
512 data points were sufficient. Due to the use of a Fast Fourier Transform subroutine, 
the number of data points must be a power of 2 (128, 256, 512, 1024, etc.). If the data 
file to be used does not have as many data points as were chosen, the program will au¬ 
tomatically pad the remaining points with zeros. The program will interpret this as an 
input signal value of zero. 

The second limiting decision which must be made is what time step size (A;) to use. 
Initially the program reads the data and calculates the average step size. From this av¬ 
erage step size the maximum time, maximum frequency, and frequency resolution are 
calculated. These quantities are related a<- follows: 


nuime — dp x At 

(16) 

A/= —-—— 

4 x miime 

(17) 

mfreq = 2 x dp x A/ 

(IS) 


where: 

miime = maximum time 

mfreq — maximum frequency 

dp = the number of data points 

A : = time step or time resolution 

A f = frequency step or frequency resolution 

In machinery monitoring applications usually two frequency spans are used, a 
broadband and a narrow band measurement. These displays may be obtained by either 
changing the filters before digitizing the data or by varying the time step size in the 
computer program. A third possibility, which has not been implemented but could be, 
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is the inclusion of a digital filter in the computer program. By selecting the smallest 
possible At the largest frequency span is obtained. The computer program will narrow 
the frequency span and focus on the lower frequencies by using a larger At. The varying 
size of the At is achieved by using every other data point or larger multiples of data 
points. Obviously the maximum frequency possible is determined by the At at 
digitization. Sometimes the At after digitization is not constant and the program allows 
for this by calculating an average At and using this average value as the actual At . 

C. WIGFUN1 

As mentioned earlier, W1GFUN1 is a time domain preprocessing program. It reads 
in the raw time data and outputs an analytic data file. Plotting subroutines are included 
in order to output various curves throughout the process. Before the analytic signal is 
calculated there are several options which may be exercised. The implications and ef¬ 
fects of each are discussed below. 

1. Effect of a Mean Value 

Figure 9 on page 19 and Figure 10 on page 20 demonstrate the effect of a mean 
value in the time domain. Figure 9 is a sine wave whose frequencies have been slowly 
varied and which has a mean value. As can be seen, the mean value in the time domain 
appears as a DC (0 Hz) component in the frequency domain. Figure 10 shows the same 
linear chirp which has had the mean value removed. As expected, the DC component 
has been eliminated. Figure 11 on page 21 is the input time signal for the linear chirp 
shown in Figure 9 and Figure 10. 
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Pseudo Wigner-Ville Distribution 
Linear Chirp 



27-MAR-1990 08:05:56.40 
512 data points 
Reduced to 256 x 128 
Smoothed 10 x 10 


Mean value removed 
Time amplified by 0.0 
Hamming window time 


Figure 10. Linear Chirp with Mean Value Removed from the Time Domain 
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2. Effect of Amplification of Time Signal Amplitude 

A subroutine has been included which amplifies the amplitude of the input sig¬ 
nal in the time domain. When this is used it definitely alters the energy distribution 
representation. The only other effect which this has on the Pseudo Wigner-Ville Dis¬ 
tribution is that it serves as a constant multiplier, vary ing the amplitudes of the distrib¬ 
ution. Figure 12 on page 23 and Figure 13 on page 24 depict a 400 Hz sine wave in 
noise which demonstrates this effect. Figure 14 on page 25 is the input time signal 
which produced Figure 12 and Figure 13. 





Figure 12. No Time Signal Amplitude 
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3. Effect of Time Signal Windows 

A subroutine has been included which will window the time domain signal. The 
reason for windowing the signal in the time domain is to bring the starting and finishing 
time values of the signal to zero. There are two windows available, a modified Hamming 
window and a Hanning window. A Generalized Hamming window is defined as: [Ref. 
25) 


window(t) = 


a + (l -ft)cos(-y-) 


0 


~ T < i < T 


otherwise 


(19) 


A Hanning window is defined as: [Ref. 26] 


window(t) * 



1 - cos 



0<t<T 


0 


otherwise 


( 20 ) 


The modified Hamming window used applied a cosine weighting function to the 
beginnning and ending of the input time record and did not vary the amplitudes inbe- 
tween. The equation used was: 

y(l-cos jjy) 0<, t < OAT 


window(t) = 


1 


o.ir</<o.9r 


1 

2 


^1 — cos 


n{T— t) \ 
OAT ) 


0.9 T<t<T 


0 


otherwise 


( 21 ) 


A modified Hamming window is the preferable window since it alters the amplitude of 
fewer data points than a Hanning window. Figure 15 on page 2S shows a Pseudo 
Wigner-Yille Distribution of a 2500 Hz sine wave with no window. Figure 16 on page 
29 is the input time domain signal. Figure 17 on page 30 shows a Pseudo Wigner-Yille 
Distribution of a 2500 Hz sine wave with a Hanning window applied in the time domain. 
Figure IS on page 31 is the modified time domain signal. Figure 19 on page 32 shows 
a Pseudo Wigner-Yille Distribution of a 2500 Hz sine wave with the modified Hamming 
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window applied in the time domain, Figure 20 on page 33 is the modified time domain 
signal. In Figure 15 it can be seen that with no window the edges of the graph at the 
start and finish times do not quite reach a zero value, Figure 17 shows that the Hanning 
window drastically changes the distribution, and Figure 19 shows that the modified 
Hamming window provides a neater distribution. 
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Pseudo Wigner-Ville Distribution 
2500 Hz sin wave 



2-APR-1990 18:20:00.51 
512 data points 
Reduced to 64 x 32 
Smoothed 10 x 10 


Mean value not removed 
Time amplified by 0.0 
No window time 


Figure 15. Pseudo Wigner-Ville Distribution with No Window 
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Pseudo Wigner-Ville Distribution 
2500 Hz sin wave 
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2-APR-1990 18:31:35.24 
512 data points 
Reduced to 64 x 32 
Smoothed 10 x 10 


Mean value not removed 
Time amplified by 0.0 
Hanning window time 
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Figure 19. Pseudo Wigner-Ville Distribution with Modified Hamming Window 
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4. Effect of Making a Real Signal Analytic 

As discussed in chapter 2, an analytic signal has no negative frequency compo¬ 
nents. By eliminating these negative frequency components the effects of aliasing are 
greatly reduced. Figure 21 on page 35 shows the aliasing problem which results when 
a real signal is not made analytic, Figure 16 on page 29 is the input 2500 Hz sine wave. 
Figure 22 on page 36 shows a time signal which has been made analytic and has had a 
modified Hamming window applied. It should be noted that the window was applied 
before the signal was made analytic, thereby saving some computation time. 
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AmpH tude 


Pseudo Wigner Distribution 
2500 Hz sin wave 



Mean value not removed 
Time amplified by 0.0 
Hamming window time 


Figure 21. Pseudo Wigner Distribution 


2-APR-1990 18:22:21.07 
512 data points 
Reduced to 64 x 32 
Smoothed 10 x 10 
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Amplitude 


Analytic Time Signal 
2500 Hz sin wave 



2-Af»R-1990 12:19:51.50 Solid llne=REAL 

Mean value not removed Dotted line-IMAGINARY 


Figure 22. Analytic Time Signal 
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D. WIG FUN IB 

WIG FUN lb is simply a program for plotting the analytic time signal. It is set ofT 
by itself in Figure 2 on page 9 since it is not required to be run, but is included as an 
added convenience. 

E. WIGFUN2 

W1GFUN2 reads in the analytic time signal from a file and calculates the Wigner- 
Ville Distribution. The Wiener-Ville Distribution is set equal to the real pan resulting 
from the complex Fourier transform of the calculated auto correlation. The algorithm 
used is based on one written by Wahl and Bolton [Ref. S] and can be expressed as: 

v 

lVDF SJ {iAujAt) = ^Re{FFT[2At CORR(i )]} (22) 

> i 


where: 

Re = Real part 

CORR{i) = The complex auto-correlation 

CORR[i) = s(J + i — 1 )s'{j — i + 1) 1 < / <; A and j < i 

CORRi /') = 0 1 < / < .V and j > i 

CORR(2.\ - / + 2) = CORR\i) A < i < 2.Y 

1. Effect of Changing the Definition of the Wigner Distribution 

As shown in chapter 2. there are several different definitions possible for re¬ 
presenting the Wigner Distribution. The computer program calculates the auto corre¬ 
lation of the signal in subroutine CORR.INCLUDE and then takes the complex FFT 
of the result. In subroutine W1GH.INCLUDE the WDF was set equal to the real part 
resulting from the FFT. Other possibilities include setting the WDF equal to the imag¬ 
inary part, setting the WDF equal to the square of the real part, or setting the WDF 
equal to the square root of the sum of the squares of the real and imaginary parts. Ac¬ 
curate representations of known signals were obtained using just the real part so this 
was the algorithm used. Additional research could be conducted in order to investigate 
other options. 

F. WIGFUN3 

WIGFUN3 is a postprocessing program with three functions. It reads in the WDT, 
reduces the output to a desired size, and applies a sliding averaging window to the 
timc-frequcncy plane. 
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1. Effect of Reducing the Wigner-Ville Distribution 

If the Wigner-Ville Distribution is not reduced and one started with 512 data 
points, the resulting distributions would be 1024 by 512 points. Graphically this is just 
too many points to visually deal with. The finest resolution which is understandable is 
256 by 12S points. The computer program allows for the following three size reductions; 
64 by 32. 12S by 64, and 256 by 128 points. (See Figure 23 on page 39, Figure 24 on 
page 40. and Figure 25 on page 41.) 
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10.2 


Pseudo Wigner-Ville Distribution 
Linear Chirp 



16-APR-1990 09:50:29.07 

512 data points Mean value removed 

Reduced to 64 x 32 Time amplified by 0.0 

Smoothed 10 x 10_Hamming window time 


Figure 23. Pseudo Wigner-Ville Distribution Reduced to 64 x 32 Points 
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Pseudo Wigner-Ville Distribution 
Linear Chirp 


16-APR-1990 11:16:55.56 
512 data points 
Reduced to 126 x 64 
Smoothed 10 x 10 


Mean value removed 
Time amplified by 0.0 
Hamming window time 



Figure 24. Pseudo Wigner-Ville Distribution Reduced to 128 x 64 Points 
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10.2 


Pseudo Wigner-Ville Distribution 
Linear Chirp 



16-APR-1990 12:17:40.91 
512 data points 
Reduced to 256 x 128 
Smoothed 10 x 10 


Mean value removed 
Time amplified by 0.0 
Hamming window time 


Figure 25. 


Pseudo Wigner-Ville Distribution Reduced to 256 x 128 Points 
















2 . Effect of Smoothing the Wigner-Ville Distribution 

As discussed in chapter 2, smoothing the Wigner-Ville Distribution deempha- 
sizes the terms which arise from mathematical calculations and emphasizes the charac¬ 
teristic events. The smoothing, or averaging, process also makes the distribution much 
more understandable (see Figure 26 on page 43 and Figure 27 on page 44). It is im¬ 
portant to note that the smoothing process uses all of the data points which were cal¬ 
culated for the Wigner-Ville Distribution but the size reduction routine only plots the 
size desired. This averaging also effects the energy representation. Figure 28 on page 
45 shows the different views and Figure 29 on page 46 shows the detailed contours of 
the reduced, but not smoothed, distribution. The sliding averaging window is based on 
an algorithm used by Wahl and Bolton (Ref. 8 ]. At a given point on the time-frequency 
plane the weighted window, hg, is expressed as: 

i __d_ jL 

hg(ij) = - - -e 200 200 (23) 

40U-*AcjA/ 

where: 

- 10 < i < 10 

- 10<j< 10 

The weighted window is applied at each of the reduced data points which are being 
plotted. 
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16-APR-1990 12:38:27.17 
512 data points 
Reduced to 256 x 128 


Mean value removed 
Time amplified by 0 0 
Hamming window time 


Figure 26. Linear Chirp with No Smoothing 
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16-APR-1990 12:24:05.01 
512 data points 
Reduced to 256 x 128 


Frequency (Hz) 


Mean value removed 
Time amplified by 0.0 
Hamming window time 


Figure 28. Different Views of Linear Chirp with No Smoothing 
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16-APR-1990 12:40:27.31 

512 data points Mean value removed 

Reduced to 256 x 128 Time amplified by 0.0 

_ _Hamming window time 


Figure 29. Detailed Contours of Linear Chirp with No Smoothing 
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G. W1GFUN4A, W1GFUN4B, AND W1GFUN4C 

The W1GFUN4 programs plot the resulting distributions. Note that the time axis 
is changed from seconds to milliseconds. \YlGFUN4a plots the 64 by 32 point distrib¬ 
utions, WlGFUN4b plots the 12S by 64 point distributions, and \\TGFUN4c plots the 
256 by 128 point distributions. 

H. OTHER CONSIDERATIONS 

1. Effect of Noise 

Figure 12 on page 23 and Figure 13 on page 24 bring up an important topic, 
which is, what is the effect of noise on the Pseudo Wigner-Ville Distribution? 
Figure 30 on page 4S is the Pseudo Wigner-Ville Distribution of two sine waves, 100 
Hz and 500 Hz, added together with minimal noise present. Figure 31 on page 49 is the 
input time signal. Figure 32 on page 50 is the Pseudo Wigner-Ville Distribution of the 
same two sine waves in a little noise, Figure 33 on page 51 is the input time signal. 
Figure 34 on page 52 is the Pseudo Wigner-Ville Distribution of the same two sine 
waves in more noise. Figure 35 on page 53 is the input time signal. By comparing these 
distributions it can be noted that as the level of noise is increased the amplitudes of the 
distributions van but the frequency content remains consistent. Each of the input time 
signal figures are made up of an upper plot from W1GFUN1 (512 data points) and a 
lower plot from a Dynamic Signal Analyzer (1024 data points) recorded at digitization. 
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Mean value removed 
Time amplified by 0.0 
Hamming window time 


11-APR-1990 13:45:01.97 
512 data points 
Reduced to 126 x 64 
Smoothed 10 x 10 


30. Pseudo Wigner-Ville Distribution with Minimal Noise 
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Pseudo Wigner-Ville Distribution 
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11-APR-1990 14:30:35.97 
512 data points 
Reduced to 128 x 64 
Smoothed 10 x 10 


Mean value removed 
Time amplified by 0.0 
Hamming window time 
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Pseudo Wigner-Ville Distribution 
100 Hz and 500 Hz sin waves in nois 



11-APR-1990 14:36:00.50 

512 data points Mean value removed 

Reduced to 128 x 64 Time amplified by 0.0 

Smoothed 10 x 10 _Hamming window time 


Figure 34. Pseudo Wigner-Ville Distribution with More Noise 
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Figure 35. Input Time Signal «ith More Noise 
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2. Distribution Interpretation 

The interpretation of the Pseudo Wigner-Ville Distribution can be very confus¬ 
ing. Before an error is suspected or if doubt is raised as to the accuracy of the programs, 
verify the input time domain signal. The Pseudo Wigner-Ville Distribution will accu¬ 
rately portray the input, however, the distribution which results may not be expected. 
Figure 36 on page 55 is the distribution from the input time signal, Figure 37 on page 
56. Figure 38 on page 57 is the distribution from the input time signal, Figure 39 on 
page 5S. From looking at the input time domain signals, the same Pseudo Wigner-Ville 
Distribution could be expected, however, they are not. The peaks in Figure 38 result 
from the extra data before and after the chirp in the time signal. 
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Pseudo Wigner-Ville Distribution 
Linear Chirp 



16-APR-1990 12:17:40.91 
512 data points 
Reduced to 256 x 128 
Smoothed 10 x 10 


Mean value removed 
Time amplified by 0.0 
Hamming window time 


Figure 36. Pseudo Wigner-Ville Distribution 


55 






Amplitude 



56 
















































































































Raw Time Signal 
Linear Chirp - A 














































































3- Swept Sine Wave Example 

Figure 40 on page 60 through Figure 45 on page 65 were produced from a sine 
wave whose frequencies were varied with a sine wave. They demonstrate the effect of 
varying the Ai and the ability to portray time dependent frequencies. Note the presence 
of the pronounced peak values located where the sweep of the frequencies are changed 
from positive to negative and negative to positive. 




mpHtude 



Figure 40. Swept Sine Wave, Frequency Span 0-1024 Hz 
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Pseudo WIgner-Vllle Distribution 
Swept Sine Wave 
Contours from zmax/6 to zmax 



14-APR-1990 15:15:01.92 

512 data points Mean value removed 

Reduced to 256 x 128 Time amplified by 0.0 

Smoothed 10 x 10__Hamming window time 


Figure 41. Detailed Contour Plot Swept Sine Wave, Frequency Span 0-1024 Hz 
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Pseudo Wigner-Ville Distribution 
Swept Sine Wave 



14-APR-1990 17:28:26.80 

512 data points Mean value removed 

Reduced to 256 x 128 Time amplified by 0.0 

Smoothed 10 x 10 __Hamming window time 


Figure 43. Swept Sine Wave, Frequency Span 0-512 Hz 
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Reduced to 256 x 128 
Smoothed 10 x 10 


Mean value removed 
Time amplified by 0 0 
Hamming window lime 


Figure 44. Detailed Contour Plot S^ept Sine Wave, Frequency Span 0-512 Hz 
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IV. PRACTICAL EXAMPLES USING PSEUDO W1GNER-V1LLE 

DISTRIBUTION 


The figures in this chapter represent data taken from an electrically driven single- 
stage centrifugal pump. The pump speed measurements are a result of data input from 
a proximity probe measuring the proximity to a gear directly coupled to the shaft. 

Figure 46 on page 67 through Figure 51 on page 72 show the pump operating at a 
steady state condition with flow rates of 10 gallons per minute and 100 gallons per 
minute. Note that since the flow rate is constant, the speed of the pump is constant for 
all time. The speed of the pump for high and low flow rates is not very different, but it 
can be seen in the detailed contour plots that for 10 GPM the speed is about 52 Hz and 
for 100 GPM the speed is about 4S Hz. 

Figure 52 through Figure 63 show the pump in transient operations. Figure 52 on 
page 73 through Figure 55 on page 76 show the pump speeding up and then steadying 
out at a constant speed. Figure 56 on page 77 through Figure 59 on page 80 show the 
pump slowing down from a constant speed. Figure 60 on page 81 through Figure 63 
on page 84 show the pump operating at a constant speed, slowing down after being 
turned off, quickly speeding up when turned back on. and steadying out again. 
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Pseudo Wigner-Ville Distribution 
Pump Speed - Steady State 10 GPM 



14-APR-1990 18:49:30.16 

512 data points Mean value removed 

Reduced to 256 x 128 Time amplified by 5.0 

Smoothed 10 x 10_Hamming window time 


Figure 46. Pump Speed Steady State 10 GPM 
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Figure 47. Detailed Contours of Steady State 10 GPM 
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Figure 48. Input Time Signal Steady State 10 GPM 
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Pseudo Wigner-Ville Distribution 
Pump - Steady State 100 GPM 



14-APR-1990 16:58:07.94 

512 data points Mean value removed 

Reduced to 256 x 128 Time amplified by 5.0 

Smoothed 10 x 10_Hamming window time 


Figure 49. Pump Speed Steady State 100 GPM 


70 






Pseudo Wlgner-Ville Distribution 
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14-APR-1990 17:00:07.08 
512 data points 
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Smoothed 10 x 10 


Mean value removed 
Time amplified by 5.0 
Hamming window time 


Figure 50. Detailed Contours of Steady State 100 GPM 
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Figure 51. Input Time Signal Steady State 100 GPM 
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Pseudo Wigner-Ville Distribution 
Pump Speed - Transient 
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14-APR-1990 17:16:26.47 
512 data points 
Reduced to 256 x 128 
Smoothed 10 x 10 


Mean value removed 
Time amplified by 5.0 
Hamming window time 


Figure 52. Pump Transient Speed lip 
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14-APR-1990 17:25:02.77 

512 data points Mean value removed 

Reduced to 256 x 128 Time amplified by 5.0 

Smoothed 10 x 10_Hamming window time 


Figure 53. Detailed Contours of Speed Up 
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14-APR-1990 17:02:50.42 

512 data points Mean value removed 

Reduced to 256 x 128 Time amplified by 5.0 

Smoothed 10 x 10_Hamming window time 


Figure 54. Multiple Views of Speed t’p 
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Pseudo Wigner-Ville Distribution 
Pump Speed - Transient 


14-APR-1990 18:42:09.83 
512 data points 
Reduced to 256 x 128 
Smoothed 10 x 10 


Mean value removed 
Time amplified by 5 0 
Hamming window time 


Figure 56. Pump Transient Coast Down 
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Figure 57. Detailed Contours of Coast Down 
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Figure 58. Multiple Views of Coast Down 
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Pseudo Wigner-Ville Distribution 
Pump Speed - Transient 



14-APR-1990 17:10:26.30 
512 data points 
Reduced to 256 x 128 
Smoothed 10 x 10 


Mean value removed 
Time amplified by 5.0 
Hamming window time 


Figure 60. Pump Transient Coast Down and Speed Up 







Figure 61. Detailed Contours or Coast Down and Speed Up 
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Figure 62. Multiple Views of Coast Down and Speed Up 
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V. CONCLUSIONS 


The Pseudo Wigner-Ville Distribution has been applied to analyzing data for use in 
machinery condition monitoring and diagnostics. From the research conducted it will 
be a valuable asset for analyzing transient machinery. The following conclusions can 
be drawn: 

• The Pseudo Wigner-Ville Distribution is ideally suited for portraying non¬ 
stationary time signals. 

• Time varying frequency components, as well as stationary frequency components, 
are evident in a single output. 

• The use of an analytic signal in calculating the Pseudo Wigner-Ville Distribution 
eliminates aliased frequency components. 

• The postprocessing smoothing process deemphasizes cross terms and provides a 
distribution directly related to the input signal. 

• The amplitude of the Pseudo Wigner-Ville Distribution varies with noise and with 
changing frequencies over time. 

• A mean value of an input time domain signal is a DC (0 Hz) component in the 
Pseudo Wigner-Ville Distribution. 

• A modified Hamming window applied to an input time domain signal starts and 
ends the signal at zero amplitude and minimizes the distortion of the Pseudo 
Wigner-Ville Distribution. 

• A time domain signal amplitude amplification, or gain, may be necessary to mini¬ 
mize roundoff error in calculating the Pseudo Wigner-Ville Distribution. 
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VI. RECOMMENDATIONS FOR FUTURE RESEARCH 


The computer program works well and the Pseudo Wigner-Ville Distribution will 
be a valuable asset to assist in machinery condition monitoring. Improvements to the 
computer program could include: 

• Increasing the number of data points being processed 

• Including a digital filter 

• Optimizing the routines used in order to increase the speed of calculations 

• Making the Pseudo Wigner-Ville Distribution an accurate signal energy represen¬ 
tation 

Areas where additional research is needed include: 

• Investigating the effect of using just the real part of the Wigner-Ville Distribution 
versus the real and imaginary parts or combinations thereof 

• Investigating other sliding averaging windows in the smoothing routine 

• Investigating the cause of the pronounced peak amplitudes evident in the swept 
sine wave examples 

Recommendations for application to machinery condition monitoring include: 

• Using an analog filter on all input data to ensure unwanted frequencies will not 
alias the frequencies of interest 
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APPENDIX A. COMPUTER CODE 


A. SUMMARY 
c 

c [rossano.thesis]symbols, include 

c 

c This is a description of the programs, subroutines, and symbols 

c used to calculate the Pseudo Wigner-Ville Distribution 

c 
c 

c FORTRAN PROGRAMS 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c SUBROUTINES 


c 


c 

AMPLIFY 


c 

Amplifies the time signal 


c 

ANALYTIC 


c 

Converts a real signal into an analytic one 

c 

CORR 


c 

Calculates the correlation 


c 

DATAIN 


c 

Reads the time signal data into 

an array 

c 

DATAOUT 


c 

Writes the WDF array to a file 


c 

DATA0UT2 


c 

Writes the RWDF and PWDF arrays 

to files 

c 

DTCALC 


c 

Calculates delta t from an input 

file or array 

c 

FFT 



c Calculates the Fast Fourier Transform 


WIGFUN1 

Reads in the real time signal, manipulates it, and 
writes the re c ulting analytic signal into a file 

WIGFUNlb 

Plots the analytic time signal 

WIGFUN2 

Reads in the analytic time signal, calculates the 
Wigner Distribution, and writes it to a file 

WIGFUN3 

Reads in the raw Wigner Distribution, manipulates 
it by reducing and smoothing, and writes the results 
to files 

WIGFUN4a 

Plots the 3 dimensional distributions which were 
reduced to (64,32) 

WIGFUN4b 

Plots the 3 dimensional distributions which were 
reduced to (128,64) 

WIGFUN4c 

Plots the 3 dimensional distributions which were 
reduced to (256,128) 
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a 


C 

c 

c 

c 

c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


HAMMING 

Applies a hamming window to a signal 

HANNING 

Applies a harming window to a signal 

MAXAMP 

Finds the maximum amplitude of a signal 

MEAN 

Calculates the mean value of a signal 

MEANR 

Removes the mean value from a signal 

MINAMP 

Finds the minimum amplitude of a signal 

PL0T2D 

Plots one 2 dimensional curve 

PL0T2D2 

Plots two 2 dimensional curves on one plot 

PL0T3D32 

Plots one 3 dimensional graph of rwdf(64,32) 

PL0T3D64 

Plots one 3 dimensional graph of rwdf(128,64) 

PL0T3D128 

Plots one 3 dimensional graph of rwdf(256,128) 
PLOT3DSPLIT32 

Plots four 3 dimensional graphs of rwdf(64,32) 
in a split screen format each with a different 
viewing perspective 
PL0T3DSPLIT64 

Plots four 3 dimensional graphs of rwdf(128,64) 
in a split screen format each with a different 
viewing perspective 
PL0T3DSPLIT128 

Plots four 3 dimensional graphs of rwdf(256,128) 
in a split screen format each with a different 
viewing perspective 
PL0TC0N32 

Plots several levels on a single contour plot for 
rwdf(64,32) 

PL0TC0N64 

Plots several levels on a single contour plot for 
rwdf(128,64) 

PL0TC0N128 

Plots several levels on a single contour plot for 
rwdf(256,128) 

PSEUDO 

Smooths the WDF along both the time and frequency axes 

RANGE32 

Finds the maximum and minimum amplitudes of 
rwdf(64,32) 

RANGE64 

Finds the maximum and minimum amplitudes of 
rwdf(128,64) 

RANGE128 

Finds the maximum and minimum amplitudes of 
rwdf(256,128) 

REDUCE 

Reduces the data in WDF to RWDF 
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SIZE 


TIMESIG 

WIGH 


WINDOW 


Changes array sizes for time plotting routines 
Modifies and plots the time signal 
Calculates the WDF 

Calls the available windowing functions 


SYMBOLS 


ain(j)=signal amplitude array 

amp=input data amplitude 

ampl=amplification applied to a signal 

amplif=amplification information for plot 

anq=answer to analytic question 

atvq=answer to various questions 

avgdelt=average dt found in input time array 

c(j)=complex correlation array 

ci(j)=imaginary correlation data array 

coef=coefficient vsed in mathematical calculations 

const=constant used in mathematical calculations 

cr(j)=real correlation data array 

datetime=date and time obtained from the computer 

delt(j)=array of actual dt values 

delta=a time increment used in calculating hamming window 

df=frequency step 

dimf=dimension size for frequency 

dimt=dimension size for time 

dp=number of data points 

dp2=dp"2 

dpend=modification of dp in input array gathering to ensure 
that dp number of points will be collected 
dpnum=number of data point information for plot 
dt.=time step 

dum-complex variable used in mathematical calculations 

dum2=complex variable used in mathematical calculations 

dum3=complex variable used in mathematical calculations 

hg(i,j)=smoothing multiplier 

i=a counter 

ii=a counter 

inname=input file name 

inv=indicator for fft or inverse fft 

iplot_val=answer whether graph output on screen or hardcopy 

j=a counter 

jj=a counter 

k=a counter 

kk=a counter 

L=a counter 

LL=a counter 

label=instructions for renaming stdOOOOl. dat 
lfile=plotting file name 
meani=mean information for plot 
meanv=mean value of a signal 
mesh=mesh size for plotting 
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c mfreq=max frequency (Hz) 

c mm=a variable used for reducing the number of time points in 

c pwdf and rwdf plots 

c mnq=answer for mean question 

c mt=number of points smoothed, time 

c mt2=mt*2 

c mtime=max time (sec) NOTE: converted to (msec) for pwdf and 

c rwdf plots 

c n=a counter 

c nf=number of points smoothed, frequency 

c nf2=nf*2 

c nn=a variable used for reducing the number of frequency points 

c in pwdf and rwdf plots 

c outname=output file name 

c pi=mathematical quantity, pi 

c pi2=2 ,v pi 

c pwdf(i,j)=pseudo wigner distribution 

c rcnum=variable for which correlation to plot 

c rdp=reduced number of data points 

c rdp2=rdp*2 

c rnq=answer for rename std00001.dat question 

c rval=reduction value, indicates number of data points 

c reduced to 

c rwdf(i,j)=reduced wigner distribution 

c s(j)=complex time array 

c si(j)=imaginary time data array 

c sizedp=maximum size available, set by dimension of arrays 

c sr(j)=real time data array 

c ss=step size used to modify dt thereby changing max freq 

c startime=a holder for datetime 

c status=variable used with datetime 

c sum=sum of values 

c sumb=sum of values 

c sval=sine of val 

c t(j)=data array 

c tcode=time window code 

c tim=input data time 

c tin(j)=signal time array 

c title=title of a plot 

c titlehold=holds 3d plot title name while plotting correlation 

c val=value used for mathematical calculations 

c wdf(i,j)=raw wigner distribution 

c windw=value of window function at a point 

c wndwcode=type of window code 

c x(j)=data array 

c xmax=raaximum amplitude 

c xmin=minimum amplitude 

c y(j)=data array 

c ymax=maximum amplitude 

c ymin=miniraum amplitude 

c zhold=holder for calculating zmax or zmin 

c zmax=maximum amplitude 

c zmin=minimum amplitude 

c 

c NOTE: The plotting routines use CA-DISSPLA Version 11. 0 
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c 

c 


written by Computer Associates International, Inc. 


B. FORTRAN PROGRAMS 

c 

c [rossano. thesis]wigfunl. for 

c 

c Graham W. Rossano 

c Prof. J. Hamilton 

c Prof. Y. S. Shin 

c 

c This program is the first of four that calculate and plot the 

c Wigner Distribution Function and variations of a signal 

c 

c KIGFUN1 reads in the real time signal, manipulates it, and writes 

c the resulting analytic signal into a file 

c 

c All of the subroutines have been broken off for ease of modifying 

c and printing. Upon compiling, they are all included, 

c 

c This program uses an include statement to facilitate changing 

c the size of the plotting arrays, 

c (see [rossano.thesis]size.include) 

c 

c A description of the symbols used in this program and its 

c subroutines may be found in [rossano.thesis]symbols, include 

c 


complex s(1025) 

real ain(1025),mtime,mfreq,avgdelt,dt,sr(1025),si(1025), 
&tin(1025),df,ampl 
character*23 datetime,tcode 
character*25 inname,Ifile 
character*50 title,outname 
character*80 label 
integer*4 status 

integer dp,dimt,dimf,rnq,dp2,j,mnq 
include 'size.include' 

c Start time of program 

status = lib$date_time (datetime) 
print* 

type *, datetime 

c Filename assignments 
print* 

print*,' Program to calculate and plot' 
print*,' Wigner Distribution Function' 
print* 

print*,' WIGFUN1' 
print* 

print*,' Enter signal input filename' 

read(5,904) inname 

print* 

print*,' Enter plotting output filename' 
read(5,901) outname 
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print* 


c Calculations 
dimt=1025 

call dtcalc (dimt,1,inname.tin,avgdeIt) 
print*,' Average Delta t = 
write(5,906) avgdelt 
print* 

print*,' Input the number of data points to evaluate' 

print*,' (128, 256, 512, or 1024)' 

print* 

print*,’ This size depends on the size of the arrays set in:' 

print*,' [rossano.thesis]size.include' 

print* 

print*,' If you change it you will need to recompile and lnkdss' 
100 read(5,*) dp 

print* 

if (dp.ne. sizedp) then 
print*,' sizedp=' 
print*,sizedp 
print*,' Reenter dp' 
go to 100 
endif 

if (dp.eq.128) then 
dimt=130 
dimf=260 
endif 

if (dp.eq.256) then 
dimt=260 
dimf=520 
endif 

if (dp. eq. 512) then 
dimt=520 
dimf=1040 
endif 

if (dp.eq.1024) then 
dimt=1025 
dimf=2050 
endif 

call datain (d mt,inname,avgdeIt,dp,ain,dt) 

tcode='No window time' 

wndwcode=0 

print*,' Input signal has been put into an array' 
print* 

call timesig (dimt,ain,dp,dt,outname,tcode, 

&title,s,sr,si,mnq,amp1) 
open (unit=7,file='wigl. log',status='new') 
open (unit=8,file='wigl. dat',status='new') 
write(7,908) dimt,dimf,dp,dt 
write(7,905) outname.tcode,mnq 
write(7,909) title,ampl 
do 200 j=l,dp 
write(8,*) s(j) 

200 continue 

close (unit=7) 
close (unit=8) 
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print*,' Do you want to rename std00001.dat?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) rnq 
print* 

if (rnq. eq. 1) then 

print*. Enter Laser plot file name' 

read(5,904) lfile 

LABEL(1:20)='RENAME STD00001.DAT ’ 

LABEL(21:45)=lfile 
CALL LIB$SPAWN(LABEL) 
print* 
end if 

print*,' Analytic data is in wigl.dat' 
print*,' Plotting info is in wigl.log' 
print* 

print*,' Now run WIGFUN2' 
print* 

c Format statements 

901 format(a50) 

904 format(a25) 

905 format(2x,a50,2x,a23,2x,i8) 

906 format(f16.8) 

908 format(2x,i8,2x,i8,2x,i8,2x,el6. 8) 

909 format(2x,a50,2x,fl6.8) 
end 

c SUBROUTINES 

include 'AMPLIFY. INCLUDE' 
include 'ANALYTIC. INCLUDE' 
include ’DATAIN. INCLUDE’ 
include 'DTCALC. INCLUDE' 
include 'HAMMING.INCLUDE' 
include 'HANNING.INCLUDE' 
include 'MAXAMP. INCLUDE' 
include 'MEAN. INCLUDE’ 
include 'MEANR. INCLUDE' 
include 'MINAMP.INCLUDE' 
include 'PL0T2D. INCLUDE' 
include 'PL0T2D2. INCLUDE' 
include 'TIMESIG. INCLUDE' 
include 'WINDOW. INCLUDE* 
c 

c [rossano. thesis]wigfunlb.for 

c 

c Graham W. Rossano 

c Prof. J. Hamilton 

c Prof. Y. S. Shin 

c 

c This program plots the analytic time signal 

c 

c All of the subroutines have been broken off for ease of modifying 

c and printing. Upon compiling, they are all included, 

c 

c A description of the symbols used in this program and its 




c subroutines may be found in [rossano.thesis]symbols.include 

c 


complex s(550) 

real mtime,mfreq,dt,sr(550),si(550),df,ampl 
character*23 datetime,tcode 
character*25 lfile 

character*50 title,outname,titlehold 
character*80 label 
integer*4 status 

integer dp,dimt,dimf,rnq,dp2,sizedp,mnq 

c Start time of program 

status = lib$date_time (datetime) 
print* 

type *, datetime 
sizedp=512 

c Get data from WIGFUN1 
print* 

print*,' Program to plot the analytic time signal' 
print* 

print*,’ WIGFUNlb' 
print* 

open (unit=7,file='wigl. log',status=’old') 
open (unit=8,file='wigl. dat',status='old') 
rewind 7 
rewind 8 

read(7,908) dimt,dimf,dp,dt 
read(7,905) outname,tcode,mnq 
read(7,909) title,ampl 
do 200 j=l,dp 
read(8,*) s(j) 

200 continue 

close (unit=7) 
close (unit=8) 
print*,' Data is loaded' 
print* 

c Calculations 

if (dp. gt.sizedp) then 

print*,' max number of points allowed is' 

print*,sizedp 

print*,' dp= 

print*,dp 

print* 

print*,' Check your dimensions' 
stop 
endif 
dp2=dp*2 
mtime=dp*dt 
df=l. /(4.*mtime) 
mfreq=2.*dp*df 
do 100 j=l,dp 
sr(j)=real(s(j)) 
si(j)=aimag(s(j)) 
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100 continue 
print* 

print*,' Do you want to plot the analytic time signal?’ 
print*,' (1 for yes, 2 for no)' 

read(5,*) ptaq 
if (ptaq.eq.1) then 
titlehold=title 
title='Analytic Time Signal' 
call plot2d2 (sr,si,dt,dp,title,outname,mnq) 
title=titlehold 
endif 

print*,' Do you want to rename std00001.dat?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) rnq 
print* 

if (rnq. ecj. 1) then 

print*, Enter Laser plot file name' 

read(5,904) lfile 

LABEL(1: 20)='RENAME STD00001.DAT ’ 

LABEL(21:45)=lfile 
CALL LIB$SPAWN(LABEL) 
endif 

print*,' WIGFUNlb. FOR is complete' 

c Format statements 

904 format(a25) 

905 format(2x,a50,2x,a23,2x,i8) 

908 format(2x,i8,2x,i8,2x,i8,2x,el6. 8) 

909 format(2x,a50,2x,fl6.8) 
end 

c SUBROUTINES 

include 'MAXAMP. INCLUDE' 
include 'MINAMP.INCLUDE' 
include 'PL0T2D2. INCLUDE' 
c 

c [ rossano. thesis]wigfun2.for 

c 

c Grahar: W. Rossano 

c Prof. J. Hamilton 

c Prof. Y. S. Shin 

c 

c This program is the 2nd of four to calculate and plot the 

c Wigner Distribution Function and variations of a signal 

c 

c WIGFUN2 reads in the analytic time signal from a file, 

c calculates the Wigner Distribution Function, and writes the 

c array to a file, 

c 

c All of the subroutines have been broken off for ease of modifying 

c and printing. Upon compiling, they are all included, 

c 

c A description of the symbols used in this program and its 

c subroutines may be found in {rossano.thesis]symbols, include 

c 
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complex s(550),c(1100) 

real mtime,mfreq,dt,cr(1100),ci(1100),df, 
&wdf(1100,550),ampl 
character*23 datetime,tcode 
character*50 title,outname 
integer ,v 4 status 

integer dp,dimt,dimf,dp2,sizedp,mnq 
common /wdfc/ wdf 

c Start time of program 

status = lib$date_time (datetime) 
print* 

type *, datetime 
sizedp=512 

c Get data from WIGFUN1 
print* 

print*,' Program to calculate and plot' 
print*,' Wigner Distribution Function' 
print* 

print*,' WIGFUN2' 
print* 

open (unit=7,file='wigl. log',status='old' ) 
open (unit=8,file='wigl. dat',status='old') 
rewind 7 
rewind 8 

read(7,908) dimt,dimf,dp,dt 
read(7,905) outname,tcode,mnq 
read(7,909) title,ampl 
do 200 j=l,dp 
read(8,*) s(j) 

200 continue 

close (unit=7) 
close (unit=8) 
print*,' Data is loaded' 
print* 

c Calculations 

if (dp.gt.sizedp) then 

print*,' max number of points allowed is' 

print*,sizedp 

print*,' dp=' 

print*,dp 

print* 

print*,' Check your dimensions' 
stop 
end if 
dp2=dp*2 
mtime=dp*dt 
df=l./(4.*mtime) 
mfreq=2.*dp*df 

call wigh (dimt,dimf,s,c,df,dt,dp2,dp) 

call dataout (dp) 

print* 

print*,' Now run WIGFUN3' 
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print* 


c Format statements 
905 format(2x,a50,2x,a23,2x,i8) 

908 format(2x,i8,2x,i8,2x,i8,2x,el6.8) 

909 format(2x,a50,2x,fl6. 8) 
end 


c SUBROUTINES 

include 'CORR.INCLUDE' 
include 'DATAOUT.INCLUDE’ 
include ’FFT.INCLUDE' 
include 'WIGH.INCLUDE' 
c 

c [rossano. thesis]wigfun3. for 

c 

c Graham W. Rossano 

c Prof. J. Hamilton 

c Prof. Y. S. Shin 

c 

c This program is the 3rd of four to calculate and plot the 

c Wigner Distribution Function and variations of a signal 

c 

c WIGFUN3 reads in the raw Wigner Distribution, manipulates 

c it by reducing and smoothing, and writes the resulting 

c arrays to files 

All of the subroutines have been broken off for ease of modifying 
c and printing. Upon compiling, they are all included, 

c 

c A description of the symbols used in this program and its 

c subroutines may be found in [ rossano. thesis]symbols, include 

c 


real mtime,mfreq,dt,df, 

&wdf(1100,550),rwdf(256,128),ampl 
character*23 datetime,tcode 
character*25 inname 
character*50 title,outname 
ir.teger*4 status 

integer dp,dimt,dimf,dp2,mnq,nn,mm,sizedp,rval 
common /wdfc/ wdf 

c Start time of program 

status = lib^date_time (datetime) 
print* 

type *, datetime 

sizedp=512 

nn=l 

mm=l 

c Get data from WIGFUN2 
print* 

print*,' Program to calculate and plot' 
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print*,' Wigner Distribution Function' 
print* 

print*,' WIGFUN3' 
print* 

open (unit=7,file='wigl. log',status='old') 
rewind 7 

read(7,908) dimt,dimf,dp,dt 
read(7,905) outname,tcode,mnq 
read(7,909) title,ampl 
close (unit=7) 

print*,' Wigl.log is loaded' 
print* 


c Calculations 

if (dp.gt.sizedp) then 

print*,' max number of points allowed is' 

print*,sizedp 

print*,' dp= 

print*,dp 

print* 

print*,' Check your dimensions' 
stop 
endif 
dp2=dp*2 
mtime=dp*dt 
df=l./(4.*mtime) 
mfreq=2. *dp*df 

open (unit=13,file='wdf. unf',status='old' , 
&forir.= 'unformatted' ) 
rewind 13 

print*,' reading wdf from wdf.unf' 
print* 

do 100 j=l,dp 
do 200 i=l,dp2 

read(13) wdf(i,j) 

200 continue 
100 continue 

close (unit=13) 

print*,' wdf array is loaded' 

print* 

call reduce (dimt,dimf,dt,dp,nn,mm,rwdf,rval) 
print* 

call dataout2 (rwdf,3,dp,nn,mm) 
print* 

call pseudo (dimt.dimf,rwdf,dp,dt,nn,mm) 
print* 

call dataout2 (rwdf,2,dp,nn,mm) 
print* 

if (title.eq.'Wigner Distribution') then 
title='Pseudo Wigner Distribution' 
endif 

if (title.eq.'Wigner-Ville Distribution') then 
title='Pseudo Wigner-Ville Distribution' 
endif 

open (unit=7,file='wig3.log',status='new') 
write(7,908) dimt,dimf,dp,dt 
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write(7,905) outname,tcode,mnq 
write(7,907) title,rval,ampl 
close (unit=7) 

print*,' Updated plotting info is in wig3.log' 
print* 

if (rval.eq.1) print*,' Now run WIGFUN4a' 
if (rval.eq.2) print*,' Now run WIGFUN4b' 
if (rval.eq. 3) print*,' Now run WIGFUN4c' 
print* 


c Format statements 
905 format(2x,a50,2x,a23,2x,i8) 

907 format(2x,a50,2x,i8,2x,fl6.8) 

908 format(2x,i8,2x,i8,2x,i8,2x,el6.8) 

909 format(2x,a50,2x,fl6. 8) 
end 

c SUBROUTINES 

include 'DATA0UT2. INCLUDE' 
include 'PSEUDO. INCLUDE' 
include 'REDUCE. INCLUDE' 
c 

c [rossano.thesis]wigfun4a.for 

c 

c Graham V. Rossano 

c Prof. J. Hamilton 

c Prof. Y. S. Shin 

c 

c This program is the 4th of four to calculate and plot the 

c Wigner Distribution Function and variations of a signal 

c 

c WIGFUN4a plots the distributions which were reduced to 

c 64 x 32 

c 

c All of the subroutines have been broken off for ease of modifying 

c and printing. Upon compiling, they are all included, 

c 

c A description of the symbols used in this program and its 

c subroutines may be found in [rossano.thesis]symbols, include 

c 


real mtime,mfreq,dt,df,rwdf(64,32),ampl 
character*23 datetime,tcode 
character*25 inname,Ifile 
character*50 title,outname,titlehold 
character*80 label 
integer*4 status 

integer dp,atvq,dimt,dimf,rnq,dp2,sizedp,n,mnq, 
&i,j,rval,rdp,rdp2 

c Start time of program 

status = lib$date_time (datetime) 
print* 

type *, datetime 
sizedp=512 
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c Get data from WIGFUN3 
print* 

print*,' Program to calculate and plot' 
print*,' Wigner Distribution Function' 
print* 

print*,' WIGFUN4a reduced to 64 x 32' 
print* 

open (unit=7,file='wig3.log',status='old') 
rewind 7 

read(7,908) dimt,dimf,dp,dt 
read(7,905) outname,tcode,mnq 
read(7,907) title,rval,ampl 
close (unit=7) 

print*,' Wig3.log is loaded' 
if (rval. eq. 1) then 
rdp=32 
else 
print* 

print*,' You are using the wrong program' 
if (rval.eq.2) print*,' Run WIGFUN4b 
if (rval.eq. 3) print*,' Run WIGFUN4c' 
stop 
endif 

if (dp.gt. sizedp) then 
print* 

print*,' max number of points allowed is' 

print*,sizedp 

print*,' dp= 

print*,dp 

print* 

print*,' Check your dimensions' 
stop 
endif 

c Calculations 
dp2=dp*2 
rdp2=rdp*2 
mtime=dp*dt 
df=l./(4. *mtime) 
mfreq=2.*dp*df 

c Convert mtime to msec for plotting 
mtime=mtime*1000. 
print* 

print*,’ Do you want to plot the reduced wdf?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
if (atvq.eq. 1) then 

open (unit=3,file='rwdf.out',status='old') 
rewind 3 
do 100 j=l,rdp 
read(3,928) n 
read(3,929) 
do 200 i=l,rdp2 

read(3,926) rwdf(i,j) 


100 




200 

continue 


read(3,929) 


read(3,929) 

100 

continue 


close (unit=3) 
titlehold=title 

if (title.eq.’Pseudo Wigner Distribution') then 
title='Reduced Wigner Distribution’ 
end if 

if (title.eq.’Pseudo Wigner-Ville Distribution’) then 
title=’Reduced Wigner-Ville Distribution’ 
end if 
print* 

print*,’ rwdf is loaded’ 
print* 

print*,’ Do you want to plot the rwdf split screen view?’ 
print*,’ (1 for yes, 2 for no)’ 

read(5,*) atvq 

if (atvq.eq.1) call plot3dsplit32 (rwdf,outname,mtime, 
&mf req, dp, tcode, t it le, mnq, ampl) 
print* 

print*,’ Do you want to plot the rwdf single view?’ 
print*,’ (1 for yes, 2 for no)’ 

read(5,*) atvq 

if (atvq.eq.1) call plot3d32 (rwdf,outname, 

&mtime,mfreq,dp,tcode,tit 1e,mnq,amp1) 
print* 

print*,’ Do you want to plot the rwdf contours?’ 
print*,’ (1 for yes, 2 for no)’ 

read(5,*) atvq 

if (atvq.eq.1) call plotcon32 (rwdf,outname, 

&mtime,mf req,dp,tcode.title,mnq,amp1) 
ticle=titlehold 
endif 
print* 

print*,’ Do you want to plot the pseudo wdf?’ 
print*,’ (1 for yes, 2 for no)’ 

read(5,*) atvq 
if (atvq.eq.1) then 

open (unit=3,file=’pwdf. out’,status=’old’) 
rewind 3 
do 300 j=l,rdp 
read(3,928) n 
read(3,929) 
do 400 i=l,rdp2 

read(3,926) rwdf(i,j) 


400 

continue 


read(3,929) 


read(3,929) 

300 

continue 


close (unit=3) 
print* 

print*,’ pwdf is loaded’ 

print* 

print* 

print*,’ Do you want to plot the pwdf split screen view?’ 
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print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plot3dsplit32 (rwdf,outname,mtime, 
&mfreq,dp,tcode,title,mnq,ampl) 
print* 

print*,' Do you want to plot the pwdf single view?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plot3d32 (rwdf,outname, 

&mtime,mfreq,dp,tcode,title,mnq,ampl) 
print* 

print*,' Do you want to plot the pwdf contours?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plotcon32 (rwdf,outname, 

&mtime,mfreq,dp,tcode,title,mnq,amp1) 
endif 

print*,' Do you want to rename std00001.dat?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) rnq 
print* 

if (rnq. eq. 1) then 

print*, Enter Laser plot file name' 

read(5,904) lfile 

LABEL(1: 20)='RENAME STD00001.DAT ' 

LABEL(21:45)=lfile 
CALL LIB$SPAWN(LABEL) 
endif 
print* 

print*,' Program is complete' 
print* 


c Format statements 

904 format(a25) 

905 format(2x,a50,2x,a23,2x,i8) 

907 format(2x,a50,2x,i8,2x,fl6.8) 

908 format(2x,i8,2x,i8,2x,i8,2x,el6.8) 

909 format(2x,a50) 

926 format(2x,el6.8) 

928 format(2x,i6) 

929 format(2x) 
end 

c SUBROUTINES 


include 'PL0T3D32. INCLUDE' 
include ’PL0T3DSPLIT32. INCLUDE' 
include 'PL0TC0N32.INCLUDE' 
include 'RANGE32.INCLUDE' 
c 

c [rossano.thesis]wigfun4b. for 

c 

c Graham W. Rossano 

c Prof. J. Hamilton 

c Prof. Y. S. Shin 

c 
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This program is the 4th of four to calculate and plot the 
Wigner Distribution Function and variations of a signal 


c 
c 
c 

c WIGFUN4b plots the distributions which were reduced to 

c 128 x 64 

c 

c All of the subroutines have been broken off for ease of modifying 

c and printing. Upon compiling, they are all included, 

c 

c A description of the symbols used in this program and its 

c subroutines may be found in [rossano. thesis]symbols, include 


real mtime,mfreq,dt,df,rwdf(128,64),ampl 
character*23 datetime,tcode 
character*25 inname,lfile 
character*50 title,outname,titlehold 
character*80 label 
integer ,v 4 status 

integer dp,atvq,dimt,dimf,rnq,dp2,sizedp,n,mnq, 
&i, j,rval,rdp,rdp2 

c Start time of program 

status = lib$date_time (datetime) 
print* 

type *, datetime 
sizedp=512 

c Get data from WIGFUN3 
print* 

print*,' Program to calculate and plot' 
print*,' Wigner Distribution Function' 
print* 

print*,' WIGFUN4b reduced to 128 x 64' 
print* 

open (unit=7,file='wig3. log',status='old') 
rewind 7 

read(7,908) dimt,dimf,dp,dt 
read(7,905) outname,tcode,mnq 
read(7,907) title,rval,ampl 
close (unit=7) 

print*,' Wig3.log is loaded' 
if (rval.eq. 2) then 
rdp=64 
else 
print* 

print*,' You are using the wrong program' 
if (rval.eq.1) print*,' Run WIGFUN4a r 
if (rval.eq.3) print*,' Run WIGFUN4c' 
stop 
endif 

if (dp.gt.sizedp) then 
print* 

print*,' max number of points allowed is' 
print*,sizedp 
print*,' dp= 
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print*,dp 
print* 

print*,' Check your dimensions' 
stop 
end if 

c Calculations 
dp2=dp*2 
rdp2=rdp*2 
mtime=dp*dt 
df=l. /(4. *mtime) 
mfreq=2. *dp*df 


c Convert mtime to msec for plotting 
mtime=mtime*1000. 
print* 

print*,' Do you want to plot the reduced wdf?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
if (atvq. eq. 1) then 

open (unit=3,file='rwdf. out’,status=’old') 
rewind 3 
do 100 j=l,rdp 
read(3,928) n 
read(3,929) 
do 200 i=l,rdp2 

read(3,926) rwdf(i,j) 

200 continue 

read(3,929) 

read(3,929) 

100 continue 

close (unit=3) 
titlehold=title 

if (title.eq.'Pseudo Wigner Distribution') then 
title='Reduced Wigner Distribution' 
endif 

if (title.eq.'Pseudo Wigner-Ville Distribution') then 
title='Reduced Wigner-Ville Distribution' 
endif 
print* 

print*,' rwdf is loaded' 
print* 

print*,' Do you want to plot the rwdf split screen view?’ 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plot3dsplit64 (rwdf,outname,mtime, 
&rafreq,dp,tcode,title,mnq,ampl) 
print* 

print*,' Do you want to plot the rwdf single view?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plot3d64 (rwdf,outname, 

&ratime,mfreq,dp,tcode.title,mnq,amp1) 
print* 

print*,' Do you want to plot the rwdf contours?' 
print*,' (1 for yes, 2 for no)' 
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read(5,*) atvq 

if (atvq.eq.1) call plotcon64 (rwdf.outname, 
Scmtime,mfreq,dp,tcode,title,mnq,ampl) 
title=titlehold 
end if 
print* 

print*,' Do you want to plot the pseudo wdf?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
if (atvq.eq. 1) then 

open (unit=3,file='pwdf. out',status='old') 
rewind 3 
do 300 j=l,rdp 
read(3,928) n 
read(3,929) 
do 400 i=l,rdp2 

read(3,926) rwdf(i,j) 

400 continue 

read(3,929) 

read(3,929) 

300 continue 

close (unit=3) 
print* 

print*,' pwdf is loaded' 

print* 

print* 

print*,' Do you want to plot the pwdf split screen view?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plot3dsplit64 (rwdf.outname,mtime, 
&mfreq,dp,tcode,title,mnq,ampl) 
print* 

print*,' Do you want to plot the pwdf single view?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plot3d64 (rwdf,outname, 

Amtime,mfreq,dp,tcode.title,mnq,ampl) 
print* 

print*,' Do you want to plot the pwdf contours?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plotcon64 (rwdf.outname, 

&mt itne,mf req, dp, tcode .title,mnq, amp 1) 
endif 

print*,' Do you want to rename std00001.dat?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) rnq 
print* 

if (rnq.e^.1) then 

print*. Enter Laser plot file name' 

read(5,904) lfile 

LABEL(1: 20)=’RENAME STD00001.DAT ’ 

LABEL(21: 45)=lfile 
CALL LIB$SPAWN(LABEL) 
endif 
print* 
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print*,' Program is complete' 
print* 


c Format statements 

904 format(a25) 

905 format(2x,a50,2x,a23,2x,i8) 

907 format(2x,a50,2x,i8,2x,fl6.8) 

908 format(2x,i8,2x,i8,2x,i8,2x,el6. 8) 

909 format(2x,a50) 

926 format(2x,el6. 8) 

928 format(2x,i6) 

929 format(2x) 
end 

c SUBROUTINES 

include 'PL0T3D64. INCLUDE' 
include 'PL0T3DSPLIT64. INCLUDE* 
include 'PL0TC0N64. INCLUDE' 
include 'RANGE64. INCLUDE’ 

[rossano. thesis]wigfun4c. for 
c 

c Graham W. Rossano 

c Prof. J. Hamilton 

c Prof. Y. S. Shin 

c 

c This program is the 4th of four tc calculate and plot the 

c Wigner Distribution Function and variations of a signal 

c 

c WIGFUN4c plots the distributions which were reduced to 

c 256 x 128 

c 

c All of the subroutines have been broken off for ease of modifying 

c and printing. Upon compiling, they are all included, 

c 

c A description of the symbols used in this program and its 

c subroutines may be found in (rossano. thesis]symbols, include 

c 


real mtime,mfreq,dt,df,rwdf(256,128),ampl 
character*23 datetime,tcode 
character*25 inname,Ifile 
character*50 title,outname,titlehold 
character*80 label 
integer*4 status 

integer dp,atvq,dimt,dimf,rnq,dp2,sizedp,n,mnq, 
&i,j,rval,rdp,rdp2 

c Start time of program 

status = lib$date_time (datetime) 
print* 

type *, datetime 
sizedp=512 

c Get data from WIGFUN3 
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print* 

print*,' Program to calculate and plot 1 
print*,' Wigner Distribution Function' 
print* 

print*,' WIGFUN4c reduced to 256 x 128' 
print* 

open (unit=7,file='wig3. log',status='old') 
rewind 7 

read(7,908) dimt.dimf,dp,dt 
read(7,905) outname,tcode,mnq 
read(7,907) title,rval,ampl 
close (unit=7) 

print*,' Wig3.log is loaded' 
if (rval.eq. 3) then 
rdp=128 
else 
print* 

print*,' You are using the wrong program' 
if (rval.eq. 1) print*,' Run WIGFUN’4a r 
if (rval.eq.2) print*,' Run WIGFUN4b' 
stop 
endif 

if (dp.gt.sizedp) then 
print* 

print*,' max number of points allowed is' 

print*,sizedp 

print*,' dp= 

print*,dp 

print* 

print*,' Check your dimensions' 
stop 
endif 

c Calculations 
dp2=dp*2 
rdp2-rdp*2 
mtime=dp*dt 
df=l./(4.*mtime) 
mfreq=2. *dp*df 


c Convert mtime to msec for plotting 
mtime=mtime*1000. 
print* 

print*,' Do you want to plot the reduced wdf?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
if (atvq.eq.1) then 

open (unit=3,file='rwdf. out'.status*'old') 
rewind 3 
do 100 j=l,rdp 
read(3,928) n 
read(3,929) 
do 200 i=l,rdp2 

read(3,926) rwdf(i,j) 

200 continue 

read(3,929) 
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read(3,929) 

100 continue 

close (unit=3) 
titlehold=title 

if (title, eq. 1 Pseudo Wigner Distribution') then 
title='Reduced Wigner Distribution' 
end if 

if (title, eq. 'Pseudo Wigner-Ville Distribution') then 
title='Reduced Wigner-Ville Distribution' 
endif 
print* 

print*,' rwdf is loaded' 
print* 

print*,' Do you want to plot the rwdf split screen view?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plot3dsplitl28 (rwdf,outnarae,mtime, 
&mfreq,dp,tcode,title,mnq,ampl) 
print* 

print*,' Do you want to plot the rwdf single view?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plot3dl28 (rwdf,outname, 

&mtime,mfreq,dp,tcode.title,mnq,ampl) 
print* 

print*,' Do you want to plot the rwdf contours?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq.1) call plotconl28 (rwdf.outname, 

&mtime,mfreq,dp,tcode,title,mnq,ampl) 
title=titlehold 
endif 
print* 

print*,' Do you want to plot the pseudo wdf?' 
print*,' (1 for yes, 2 for no)’ 

read(5,*) atvq 
if (atvq.eq.1) then 

open (unit=3,file='pwdf.out',status='old') 
rewind 3 
do 300 j=l,rdp 
read(3,928) n 
read(3,929) 
do 400 i=l,rdp2 

read(3,926) rwdf(i,j) 

400 continue 

read(3,929) 
read(3,929) 

300 continue 

close (unit=3) 
print* 

print*,' pwdf is loaded' 

print* 

print* 

print*,' Do you want to plot the pwdf split screen view?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
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if (atvq.eq. 1) call plot3dsplitl28 (rwdf,outname,mtime, 
&mfreq,dp,tcode,title,mnq,ampl) 
print* 

print*,’ Do you want to plot the pwdf single view?’ 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 

if (atvq.eq. 1) call plot3dl28 (rwdf,outname, 

&mtime,mfreq,dp,tcode,title,mnq.ampl) 
print* 

print*,' Do you want to plot the pwdf contours?' 
print*,' (1 for yes, 2 for no)’ 

read(5,*) atvq 

if (atvq.eq.1) call plotconl28 (rwdf,outname, 

&mtime,mf req,dp,tcode,title,mnq,amp1) 
end if 

print*,' Do you want to rename std00001.dat?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) rnq 
print* 

if (rnq.eq.1) then 

print*, Enter Laser plot file name' 

read(5,904) lfile 

LABEL(1:20)='RENAME STD00001.DAT ' 

LABEL(21: 45)=lfile 
CALL LIB$SPAWN(LABEL) 
end if 
print* 

print*,' Program is complete' 
print* 

c Format statements 

904 format(a25) 

905 format(2x,a5C,2x,a23,2x,i8) 

907 format(2x,a50,2x,i8,2x,fl6.8) 

908 format(2x,i8,2x,i8,2x,i8,2x,el6. 8) 

909 format(2x,a50) 

926 format(2x,el6. 8) 

928 format(2x,i6) 

929 format(2x) 
end 

c SUBROUTINES 

include 'PL0T3D128. INCLUDE’ 
include 'PLOT3DSPLIT128. INCLUDE' 
include ’PLOTCON128. INCLUDE’ 
include 'RANGE128. INCLUDE’ 

C. ALPHANUMERIC LISTING OF SUBROUTINES 

c 

c [rossano.thesis]amplify, include 

c 

c This subroutine amplifies the signal ain 


subroutine amplify (dimt,ain,dp,ampl) 
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integer dp,j,dimt 
real ain(dimt),amp1 

print*, ' What value would you like to amplify the signal by?' 
print*, ' (Do not forget the decimal point)' 

read(5,901) ampl 
do 100 j=l,dp 

ain(j)=ampl*ain(j) 

100 continue 
return 

901 format(f16.8) 
end 
c 

c [rossano. thesis]analytic.include 

c 

c This subroutine converts a real signal into an analytic one 


subroutine analytic (dimt,sr,s,dp) 

integer dimt,dp,i,j 
complex s(dimt) 

real sr(dimt),dt,pi,sum,val,n,sumb,sval 

pi=4. *atan( 1. ) 

do 100 i=l,dp 
sum=0. 

do 200 j=l,dp 
sumb=0. 

if(i-j.eq.0) go to 200 

n=i-j 

val=pi*n/2. 
sval=sin(val) 
sumb=sr(j)*sval*sval/val 
200 sum=sum+sumb 

100 s(i)=cmplx(sr(i),sum) 
return 
end 
c 

c [rossano.thesis]corr. include 

c 

c This subroutine calculates the correlation 


subroutine corr (dimt,dimf,s,j,dt,c,dp) 

integer dimt.dimf,dp,dp2,i,j 
complex s(dimt),c(dimf),dum 
real coef,dt 

dp2=dp*2 
coef=2.*dt 
do 100 i=l,dp+l 
if(j•ge.i) then 
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dum=s(j-i+1) 
else 

dum=cmplx(0.,0.) 
end if 

c(i)=coef*(s(j+i-l)*conjg(dum)) 
if( i. ne. j. and. i. ne. dp+1) then 
c(dp2-i+2)=conjg(c(i)) 
endif 

100 continue 
return 
end 
c 

c [rossano.thesis]datain. include 

c 

c This subroutine reads the time signal data into an array 

c 

c ASSUME input file in format: time amplitude (2x,el6.8,5x,el6.8) 
c ASSUME end of file indication is a last entry of 9999. ,9999. 
c Correct format may be obtained using [rossano.data]convert.for 
c 


subroutine datain (dimt,inname,avgdelt,dp,ain,dt) 
integer j,i,dp,ss,dpend,n,dimt 

real tin(2050),ain(dimt),amp,tim,mtime,mfreq,avgdelt,df,dt 
character*25 inname 

open (unit=4,file=inname,status='old') 
rewind 4 

mtime=dp*avgde1t 
df=l./(4. *mtime) 
mfreq=2. *dp*df 

print*,' Max time is' 

write(5.906) mtime 

print*, Delta f is' 

write(5j906) df 

print*. Max frequency is now' 

write(5,906) mfreq 

print* 

print*,' This can be changed by multiplying the delta t step size' 
print*,' if step size = 1 ; Max freq remains' 
write(5,906) mfreq 

print*,' if step size = 2 ; Max freq becomes df= mtime=' 

mtime=dp*2. *avgdelt 

df=l./(4.*mtime) 

mfreq=2. *dp*df 

write(5,907) mfreq,df,mtime 

print*,' If step size = 5 ; Max freq becomes df= mtime=' 

mtime=dp*5*avgdelt 

df=l./(4. *mtime) 

mfreq=2.*dp*df 

write(5,907) mfreq,df,mtime 
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df= 


mtime 


print*,' If step size = 10 ; Max freq becomes 
mtime=dp*10*avgdeIt 
df=l./(4.*mtime) 
mfreq=2.*dp*df 
vrite(5,907) mfreq,df,mtime 

print*,' If step size = 16 ; Max freq becomes df= 

mtime=dp*16*avgdelt 

df=l. /(4. *mtime) 

mfreq=2. *dp*df 

write(5,907) mfreq,df,mtime 

print* 

print*,' Input desired step size' 

read(*,*) ss 

print* 

do 111 j=l,dimt 
ain(j)=0.0 
tin(j)=0. 0 
111 continue 

i=0 

n=ss-1 
dpend=dp*ss 
do 200 j=l,dpend 
n=n+l 

read(4,902) tim,amp 

if (tim.eq.9999.) then 

if (amp. eq. 9999. ) go to 210 
endif 

if (n.eq.ss) then 
i=i+l 

ain(i)=amp 
tin(i)=tim 
n=0 
endif 

200 continue 
go to 212 

c This pads zeros if the signal is too short 

210 j=i+l 

print*,' Padding with zeros since not enough data points' 
print* 

do 211 j=j,dp 
ain(j)=0. 0 
tin(j)=j*ss*avgdelt 

211 continue 

212 close (unit=4) 
mtime=tin(dp) 
dt=mtime/(dp-l) 

call dtcalc (dimt,dp,inname,tin,avgdelt) 


mtime= 
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print*,' Average dt=' 
write(5.906) avgdelt 
print*, Final dt=' 
write(5,906) dt 
print* 

if (avgdelt.ne.dt) then 

print*,' Average dt does not equal Final dt' 
print*,' This is probably because your time record' 
print*,' did not start at zero.' 

print*,' This can be fixed by letting dt=average dt' 
print*,' Is this ok?' 

print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
if (atvq.eq.2) stop 
dt=avgdelt 
endif 
return 

c Format statements 
902 format(2x,el6.8,5x,el6.8) 

906 format(f16.8) 

907 format(20x,f16.8,fl6.8,fl6.8) 
end 

c 

c [rossano.thesis]dstaout.include 

c 

c This subroutine prints the WDF array to a file 

c 

subroutine dataout (dp) 

integer j,i,dp,dp2 
real wdf(1100,550) 
common /wdfc/ wdf 

open (unit=7,file='wdf. unf',status='new', 

&form='unformatted') 

print*,' Writing to wdf.unf’ 

dp2=dp*2 
do 100 j=l,dp 
do 200 i=l,dp2 

write(7) wdf(i,j) 

..00 continue 

100 continue 

close (unit=7) 
return 
end 
c 

c [rossano.thesis]dataout2.include 

c 

c This subroutine prints the RWDF and PWDF arrays to files 

c 


subroutine dataout2 (rwdf,n,dp,nn,mm) 
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integer j,i,n,dp,dp2,nn,mm 
real rwdf (256,128'' 

if (n. eq. 2) then 

open (unit=7,file"'pwdf. out'.status='new’) 
print*,' Writing to pwdf. out 
endif 

if (n.eq.3) then 

open (unit=7,file='rwdf. out',status='new') 
print*,' Writing to rwdf.out' 
endif 
dp2=dp*2 

do 100 j=l,dp/mm 
write(7,908) j 
write(7,909) 
do 200 i=l,dp2/nn 

write(7,906) rwdf(i,j) 


200 

continue 

W’rite( 7,909) 
write(7,909) 

100 

continue 
close (unit=7) 
return 

906 

format (2x,el6.8) 

908 

format (2x,i6) 

909 

format (2x) 
end 


c 

c [rossano.thesis]dtcalc.include 

c 

c This subroutine calculates delta t from an input file (j=l) or 

c an array (j=dp) 

c 


subroutine dtcalc (dimt,j,inname,tin,avgdelt) 
integer i,j,n,dimt 

real tin(dimt),ain(1025),sum,delt(1025).avgdelt 
character*25 inname 


n=0 

sum=0. 

if (j.ne.1) then 
n=j 

go to 300 
endif 

open (unit=4,file=inname,status='old') 
rewind 4 

do 100 i=l,1025 

read(4,904) tin(i),ain(i) 
if (tin(i).eq.9999.) then 

if (ain(i).eq.9999.) go to 200 
endif 
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n=n+l 

100 continue 
200 close (unit=4) 

300 do 400 i=l,n-l 

delt(i)=tin(i+l)-tin(i) 
sum=sum+delt(i) 

400 continue 

avgdelt=sum/(n-l) 

return 

904 format (2x,el6.8,5x,el6.8) 
end 
c 

c [rossano.thesis]fft.include 

c 

c This subroutine calculates the Fast Fourier Transform 

c (FFT if inv=0, Inverse FFT if inv=l) 

c 

subroutine fft (dimf,c,dp2,inv) 

integer dimf,inv,dp,dp2,i,j,n,val,ii,k 
complex c(dimf),dum,dum2,dum3 
real pi,const,coef 

pi=4.*atan(1. ) 
const=dp2 

val=alog(const)/alog(2. )+. 1 
dp=dp2/2 

j=l 

do 40 i=l,dp2-l 

if (i. ge. j) go to 10 
dum3=c(j) 
c(j)=c(i) 
c(i)=dum3 
10 k=dp 

20 if (k. ge. j) go to 30 

j=j'k 
k=k/2 
go to 20 
30 j=j+k 

40 continue 

do 70 n=l,val 
coef=2**n 
coef=coef/2 
dum2=cmplx(1. ,0. ) 

dum=cmplx(cos(pi/coef),-sin(pi/coef)) 
if (inv.ne.0) dum=conjg(dura) 
do 60 j=l,coef 

do 50 i=j,dp2,2*coef 
ii=i+coef 
dum3=c(ii)*dum2 
c(ii)=c(i)-dum3 
c(i)=c(i)+dum3 
50 continue 
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dum2=dum2 lV dum 

continue 

continue 

if (inv.eq.0) return 
do 80 i=l,dp2 

c(i)=c(i)/cmplx(const,0.) 
continue 
return 
end 

[ rossano. thesis]hamming, include 

This subroutine applies a hamming window to the signal ain 


subroutine hamming (dimt,ain,dp,dt) 
integer j,dp,dimt 

real ain(dimt),dt,mtime,pi,const,delta,n 

pi=4. 0*atan( 1.0) 
mtime=dp*dt 
delta=0.l*mtime 
const=pi/delta 

do 100 j=l,dp 
n=(j-l)*dt 

if (n. le.delta) ain(j)=ain(j)*(. 5*(1. -cos(const*n))) 
if (n.ge.mtime-delta) then 

ain(j)=ain(j)*(. 5*(1. -cos(const*(mtime-n)))) 
endif 

100 continue 
return 
end 
c 

c [ rossano. thesis]hanning. include 

c 

c This subroutine applies a hanning window to signal ain 

c 

subroutine hanning (dimt,ain,dp,dt) 
integer j,dp,dimt 

real ain(dimt),dt,mtime,pi2,windw,const 

pi2=2. *4. 0 ,v atan( 1. 0) 

mtime=dp*dt 

con s t=(p i 2 ,v d t) / m t ime 

do 100 j=0,dp-l 

windw=0. 5*( 1 -cos( const ,v j)) 
ain( j+l)=ain( j+l) ,v windw 
100 continue 
return 
end 
c 

c [rossano.thesis]maxamp.include 
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c 

c 

c 


This subroutine finds the maximum amplitude of array y 


subroutine maxamp (y,dp,zmax) 

integer i,dp 

real y(dp),zmax,zhold 

zmax=0.0 

do 200 i=l,dp 
zhold=y(i) 

if (zhold. gt. zmax) then 
zmax=zhold 
endif 

200 continue 

return 
end 
c 

c [rossano.thesis]mean, include 

c 

c This subroutine calculates the mean value in the signal ain 


subroutine mean (dimt,ain,dp,meanv) 

integer dp,i,dimt 
real ain(dimt),sum,meanv 

sum=0.0 
do 100 i=l,dp 
s'im=sum+ain( i) 

100 continue 

meanv=sum/dp 

return 

end 

c 

c [rossano.thesis]meanr.include 

c 

c This subroutine removes the mean value from the signal ain 


subroutine meanr (dimt,ain,dp,meanv) 

integer dp,i,dimt 
real ain(dimt),meanv 

do 100 i=l,dp 

ain(i)=ain(i)-meanv 
continue 
return 
end 

[rossano.thesis]minamp.include 
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u o 


c 

c 


This subroutine finds the minimum amplitude of array y 


subroutine minamp (y,dp,zmin) 

integer i,dp 

real y(dp),zmin,zhold 

zmin=0.0 

do 200 i=l,dp 
zhold=y(i) 

if (zhold.It.zmin) then 
zmin=zhold 
end if 

200 continue 
return 
end 

[rossano.thesis]plot2d. include 

c This subroutine uses disspla 11.0 to plot a 2 Dimension plot 

c 

c The data input array is y, the x axis assumes a constant dt spacing 


subroutine plot2d (y,dt,dp,title,outname) 

real dt,mtime,zmax,zmin 
character*23 datetime 
character*50 title.outname 
integer*4 status 
integer iplot_val,dp,mesh,j 
include 1 size.include' 

call maxamp (y,dp,zmax) 
call minamp (y,dp,zmin) 

mtime=dp*dt 
do 100 j=l,dp 
x(j)=(j-l)*dt 
100 continue 

print* 

print*, Enter limits for plotting (do not forget the decimal) 

print*,' YMIN YMAX' 

write(5,*) zmin,zmax 

read(5,*) zmin,zmax 

print* 

mesh=l 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
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print* 

if (iplot_val . eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
end if 

if (iplot_val .eq. 2) then 

print*,' Please be patient' 
print*,' This will take several minutes.' 
call LN03I 
endif 

call swissm 

call hwshd 

call chrpat (16) 

call height (.325) 

call physor (1.0,0.625) 

call area2d (6.0,8.5) 

call alnmes (.5,0.) 

call messag (title,100,3.2,9.9) 

call messag (outname,100,3. 2,9. 5) 

call reset ('alnmes') 

call blsur 

call height (0.200) 

call xname ('Time (sec)',10) 

call yname ('Amplitude',9) 

call cross 

call xintax 

call yintax 

call graf (0. ,mtime/4. ,mtime,zmin,(zmax-zmin)/4. ,zmax) 

print*,' Axes are complete' 

call reset ('hwshd') 

call curve (x,y,dp,0) 

call height (.150) 

call alnmes (0. 0,1.0) 

status = lib$date_time (datetime) 

call messag (datetime,23,-0. 8,-0. 4) 

call reset ('alnmes') 

print*,' Surface is complete' 

510 call end3gr (0) 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [rossano.thesis]plot2d2. include 

c 

c This subroutine uses disspla 11.0 to plot 2 2 Dimension curves 

c on one plot 

c 

c The data input arrays are x and y, the x axis assumes a constant dt 
c spacing 
c 

subroutine plot2d2 (x,y,dt,dp,title,outname,mnq) 
real dt,mtime,zmax,zmin,xmax,xmin,ymax,ymin 
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character*23 datetime 
character*50 title.outname,meani 
integer*4 status 
integer iplot_val,dp,mesh,j,mnq 
include 1 size.include' 


mtime=dp*dt 
mesh=l 

if (mnq. eq. 1) then 

meani='Mean value removed' 
end if 

if (mnq.eq. 2) then 

meani='Mean value not removed' 
end if 

call maxamp (x,dp,zmax) 
call minamp (x,dp,zmin) 
xmax=zmax 
xmin=zmin 

call maxamp (y,dp,zmax) 
call minamp (y,dp,zmin) 
ymax=zmax 
ymin=zmin 

if (xmax. gt. ymax) zmax=xmax 
if (xmin. It. ymin) zmin=xmin 


print* 

print*,' Enter limits for plotting (do not forget the decimal)’ 

print*,' YMIN YMAX’ 

write(5,*) zmin.zmax 

read(5,*) zmin.zmax 

print* 

do 100 j=l,dp 
t(j)=(j-l)*dt 
100 continue 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val . eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val .eq. 2) then 

print*,' Please be patient' 
print*,' This will take several minutes.' 
call LN03I 
endif 

call swissm 
call hwshd 
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call chrpat (16) 

call height (.325) 

call physor (1.0,0.625) 

call area2d (6.0,8.5) 

call alnmes (.5,0.) 

call messag (title,100,3.2,9.9) 

call messag (outname,100,3.2,9.5) 

call reset ('alnmes') 

call blsur 

call height (0.200) 

call xname ('Time (sec)',10) 

call yname ('Amplitude',9) 

call cross 

call xintax 

call yintax 

call graf (0.,mtime/4.,mtime,zmin,(zmax-zmin)/4.,zmax) 

print*,' Axes are complete' 

call reset ('hwshd') 

call curve (t,x,dp,0) 

call dot 

call curve (t,y,dp,0) 

call height (.150) 

call alnmes (0. 0,1.0) 

status = lib$date_time (datetime) 

call messag (datetime,23,-0.8,-0.1) 

call messag (meani,50,-0.8,-0.4) 

call messag ('Solid line=REAL',15,3.0,-0.1) 

call messag ('Dotted line=IMAGINARY',21,3.0,-0.4) 

call reset ('alnmes') 

print*,' Surface is complete' 

510 call end3gr (0) 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 

c 

c [ rossano. thesis]plot3d32. include 

c 

c This subroutine uses disspla 11.0 to plot a single 3 dimension 

c graph of rwdf(64,32) 


subroutine plot3d32 (rwdf,outname,mtime,mfreq,dp, 
6tt code, t it 1 e, mnq, amp 1) 

real mtime,mfreq,zmax,zmin,ampl 
character*23,datetime,tcode,dpnum,amp1if 
character*50.title,outname,meani 
integer*4 status 

integer iplot_va1,dp,mesh,rdp2,rdp,mnq 
real rwdf(64,32) 

rdp=32 

if (dp.eq.128) then 

dpnum='l28 data points' 
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endif 

if (dp. eq. 256) then 

dpnum- 256 data points' 
endif 

if (dp.eq.512) then 

dpnum- 512 data points' 
endif 

if (dp.eq.1024) then 

dpnum- 1024 data points' 
endif 

if (mnq. eq. 1) then 

meani='Mean value removed' 
endif 

if (mnq. eq. 2) then 

meani='Mean value not removed' 
endif 

amplif='Time amplified' 

if (ampl.eq.0.0) amplif='Time amplified by 0.0' 
if (ampl.eq.1.0) amplif='Time amplified by 1.0' 
if (ampl.eq.2.0) amplif='Time amplified by 2.0' 
if (ampl.eq.3.0) amplif='Time amplified by 3.0' 
if (ampl.eq.4.0) amplif='Time amplified by 4.0' 
if (ampl.eq.5.0) amplif='Time amplified by 5.0' 
call range32 (rwdf,rdp,zmax,zmin) 
print* 

print*,' Enter limits for plotting (do not forget decimal)' 

print*,' ZMIN ZMAX' 

write(5,*) zmin,zmax 

read(5,*) zmin.zmax 

print* 

mesh=l 

rdp2=rdp*2 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val .eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val . eq. 2) then 

print*,' Please be patient' 
print*,' This will take several minutes.' 
call LN03I 
endif 

call swissm 

call hwshd 

call chrpat (16) 

call height (.325) 

call physor (.75,.625) 

call area2d (7.5,9.75) 

call alnmes (.5,0.) 

call messag (title,100,7.5/2.,9.25) 
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call messag (outname,100,7.5/2.,8.75) 
call reset ('airlines') 
call blsur 

call volm3d (8. ,8. ,9. ) 

call x3name ('Frequency (Hz)',15) 

call y3name ('Time (msec)',11) 

call z3name ('Amplitude',9) 

call vuangl (-60. ,30. ,30. ) 

call xintax 

call yintax 

call zintax 

call graf3d (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3.,zmax) 
print*,' Axes are complete' 
call reset ('hwshd') 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 

call reset ('alnmes') 

call height (0.150) 

call alnmes (0. 0,1.0) 

call messag (meani,50,5.0,0.0) 

call messag (amplif,23,5. 0,-0. 2) 

call mess'® . tcode,23,5.0,-0.4) 

status = 1ib$date_time (datetime) 

call eie'^ag (datetime,23,0. 0,0. 2) 

call messag (dpnum,23,0.0,0. 0) 

call messag ('Reduced to 64 x 32' ,23,0.0,-0.2) 

call messag ('Smoothed 10 x 10',23,0.0,-0.4) 

call reset ('alnmes') 

print*,' Surface is complete' 

510 call end3gr (0) 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [rossano. thesis]plot3d64. include 

c 

c This subroutine uses disspla 11.0 to plot a single 3 dimension 

c graph of rwdf(128,64) 


subroutine plot3d64 (rwdf,outname,mtime,mfreq,dp, 
&tcode,title,mnq,ampl) 

real mtime,mfreq,zmax,zmin,ampl 
character*23,datetime,tcode,dpnum,amplif 
character*50.title,outname,meani 
integer*4 status 

integer iplot_val,dp,mesh,rdp2,rdp,mnq 
real rwdf(128,64) 

rdp=64 

if (dp.eq.128) then 

dpnum='l28 data points' 
end if 
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if (dp.e^. 256) then 

dpnum-256 data points' 
end if 

if (dp. ecj. 512) then 

dpnum-512 data points' 
endif 

if (dp.eq. 1024) then 

dpnum-1024 data points' 
endif 

if (mnq.eq. 1) then 

meani='Mean value removed' 
endif 

if (mnq. eq. 2) then 

meani='Mean value not removed' 
endif 

amplif='Time amplified' 

if (ampl.eq.0.0) amplif='Time amplified by 0. O' 
if (ampl.eq.1.0) amplif='Time amplified by 1. O' 
if (ampl.eq.2.0) amplif='Time amplified by 2.0' 
if (ampl.eq. 3. 0) amplif='Time amplified by 3.0' 
if (ampl.eq.4.0) amplif='Time amplified by 4.0' 
if (ampl.eq.5.0) amplif='Time amplified by 5.0' 
call range64 (rwdf,rdp,zmax,zmin) 
print* 

print*,' Enter limits for plotting (do not forget decimal)' 

print*,' ZMIN ZMAX' 

write(5,*) zmin,zmax 

read(5,*) zmin.zmax 

print* 

mesh=l 

rdp2=rdp*2 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val . eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val .eq. 2) then 

print*,' Please be patient' 
print*,' This will take several minutes.' 
call LN03I 
endif 

call swissm 
call hwshd 


call chrpat (16) 

call height (.325) 

call physor (.75,.625) 

call area2d (7.5,9.75) 

call alnmes (.5,0.) 

call messag (title,100,7.5/2.,9.25) 

call messag (outname,100,7.5/2.,8.75) 
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call reset ('alnmes') 
call blsur 

call volm3d (8. ,8. ,9. ) 

call x3name ('Frequency (Hz)',15) 

call y3name (’Time (msec)’,11) 

call z3name ('Amplitude',9) 

call vuangl (-60.,30.,30. ) 

call xintax 

call yintax 

call zintax 

call graf3d (0.,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,' Axes are complete' 
call reset ('hwshd') 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 

call reset ('alnmes') 

call height (0.150) 

call alnmes (0. 0,1.0) 

call messag (meani,50,5. 0,0. 0) 

call messag (amplif,23,5.0,-0.2) 

call messag (tcode,23,5.0,-0.4) 

status = lib$date_time (datetime) 

call messag (datetime,23,0.0,0.2) 

call messag (dpnum,23,0.0,0.0) 

call messag ('Reduced to 128 x 64',23,0.0,-0.2) 

call messag ('Smoothed 10 x 10',23,0. 0,-0.4) 

call reset ('alnmes') 

print*,' Surface is complete' 

510 call end3gr (0) 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [ rossano. thesis]plot3dl28.include 

c 

c This subroutine uses disspla 11.0 to plot a single 3 dimension 

c graph of rwdf(256,128) 


subroutine plot3dl28 (rwdf,outname,mtime,mfreq,dp, 
&tcode,title,mnq,ampl) 

real mtime,mfreq,zmax,zmin,ampl 
character*23,datetime,tcode,dpnum,amplif 
character*50,title,outname,meani 
integer*4 status 

integer iplot_val,dp,mesh,rdp2,rdp,mnq 
real rwdf(256,128) 

rdp=128 

if (dp.eq. 128) then 

dpnum='l28 data points' 
endif 

if (dp.eq.256) then 
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dpnum='256 data points' 
end if 

if (dp. eq.512) then 

dpnum- 512 data points' 
end if 

if (dp, eq. 1024) then 

dpnum-1024 data points' 
endif 

if (mnq. eq. 1) then 

meani='Mean value removed' 
endif 

if (mnq.eq. 2) then 

meani='Mean value not removed' 
endif 

amplif='Time amplified' 
if (ampl. eq. 0. 0) amplif= 
if (ampl. eq. 1. 0) amplif= 
if (ampl. eq. 2. 0) amplif= 
if (ampl. eq. 3. 0) amplif= 
if (ampl.eq. 4. 0) amplify 
if (ampl.eq.5. 0) amplif= 


'Time amplified by 
'Time amplified by 
'Time amplified by 
'Time amplified by 
'Time amplified by 
'Time amplified by 


call rangel28 (rwdf,rdp,zmax,zmin) 
print* 

print*,' Enter limits for plotting (do not 

print*,' ZMIN ZMAX' 

write(5,*) zmin,zmax 

read(5,*) zmin,zmax 

print* 

mesh=l 


0 . O' 
1. 0’ 
2 . 0 ’ 

3. O' 

4. 0' 

5. O' 


forget decimal)' 


rdp2=rdp*2 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val .eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val .eq. 2) then 

print*,' Please be patient' 
print*,' This will take several minutes.' 
call LN03I 
endif 

call swissm 

call hwshd 

call chrpat (16) 

call height (.325) 

call physor (.75,.625) 

call area2d (7.5,9.75) 

call alnmes (.5,0.) 

call messag (title,100,7.5/2. ,9.25) 

call messag (outname,100,7.5/2.,8.75) 

call reset ('alnmes') 
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call blsur 

call volm3d (8. ,8. ,9. ) 

call x3name ('Frequency (Hz)',15) 

call y3name ('Time (msec)',11) 

call z3name ('Amplitude',9) 

call vuangl (-60.,30.,30.) 

cal] xintax 

call yintax 

call zintax 

call graf3d (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3.,zmax) 
print*,' Axes are complete' 
call reset ('hwshd') 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 

call reset ('alnmes') 

call height (0.150) 

call alnmes (0. 0,1.0) 

call messag (meani,50,5.0,0. 0) 

call messag (amplif,23,5.0,-0.2) 

call messag (tcode,23,5.0,-0.4) 

status = lib$date_time (datetime) 

call messag (datetime,23,0.0,0.2) 

call messag (dpnum,23,0.0,0.0) 

call messag ('Reduced to 256 x 128',23,0.0,-0.2) 
call messag ('Smoothed 10 x 10*,23,0.0,-0.4) 
call reset ('alnmes') 
print*,' Surface is complete' 

510 call end3gr (0) 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [rossano.thesis]plot3dsplit32. include 

c 

c This subroutine uses disspla 11.0 to plot a series of 

c 3 dimensional graphs of rwdf(64,32) 


subroutine plot3dsplit32 (rwdf,outname,mtime,mfreq,dp,tcode, 
&title,mnq,ampl) 

real mtime,mfreq,zmax,zmin,rwdf(64,32),ampl 
character*23,datetime,tcode,dpnum,amplif 
chart cter*50,title,outname,meani 
integer*4 status 

integer iplot_val,dp,mesh,rdp2,rdp,mnq 
rdp=32 

if (dp.eq.128) then 

dpnum='l28 data points' 
endif 

if (dp.eq.256) then 

dpnum='256 data points' 
endif 
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if (dp.eq.512) then 

dpnum- 512 data points' 
endif 

if (dp.eq.1024) then 

dpnum- 1024 data points' 
endif 

if (mnq. eq. 1) then 

meani='Mean value removed' 
endif 

if (mnq.eq.2) then 

meani='Mean value not removed' 
endif 

amplif='Time amplified' 
if (ampl. eq. 0. 0) amplif= 
if (ampl. eq. 1. 0) amplif= 
if (ampl.eq. 2. 0) amplif= 
if (ampl. eq. 3. 0) amplif= 
if (ampl. eq. 4. 0) amplif= 
if (ampl.eq. 5. 0) amplif= 


Time amplified by 0.0' 
Time amplified by 1.0' 
Time amplified by 2.0' 
Time amplified by 3.0' 
Time amplified by 4.0' 


Time amplified by 5.0' 
call range32 (rwdf,rdp,zmax,zmin) 
print*,' Enter limits for plotting (do not forget decimal)' 
print*,' ZMIN ZMAX r 

write(5,*) zmin,zmax 
read(5,*) zmin,zmax 
print* 
mesh=l 
rdp2=rdp*2 
call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)’ 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val . eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val . eq. 2) then 

print*,' Please be patient' 
print*,' This will take several minutes.' 
call LN03I 
endif 

call swissm 
call hwshd 
call chrpat (16) 


c LABELS 

call height (0.200) 

call physor (0.5,0.625) 

call area2d (7.5,9.75) 

call alnmes (0.5,0.0) 

call raessag (title,100,7.5/2. ,9.75) 

call messag (outname,100,7. 5/2. ,9. 5) 

call reset ('alnmes') 

call height (0.150) 
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call alnmes (0.0,1.0) 

call messag (meani,50,5. 0,0. 0) 

call messag (amplif,23,5.0,-0.2) 

call messag (tcode,23,5.0,-0.4) 

status = lib$date_time (datetime) 

call messag (datetime,23,0. 0,0. 2) 

call messag (dpnum,23,0.0,0.0) 

call messag ('Reduced to 64 x 32',23,0.0,-0.2) 

call messag ('Smoothed 10 x 10',23,0.0,-0.4) 

call reset ('alnmes') 

call endgr (1) 

print*, 1 Labels are complete' 

c FIRST SUBPLOT 

call height (.150) 
call physor (.5,5.5) 
call area2d (3.5,4.875) 
cal 1 blsur 

call volm3d (8. ,8. ,9. ) 

call x3name ('Frequency (Hz)',15) 

call y3name ('Time (msec)',11) 

call z3name ('Amplitude',9) 

call vuangl (-50. ,0. ,30. ) 

call xintax 

call yintax 

call zintax 

call graf3d (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
ot(zmax-zmin)/3. ,zmax) 
print*,' First subplot axes are complete' 
call reset ('hwshd ; 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 
call endgr (2) 

print*,' First subplot surface is complete' 

c SECOND SUBPLOT 

call height (.150) 
call phvsor (4.25,5.5) 
call area2d (3.5,4.5) 
call blsur 

call volm3d (8. ,8. ,9. ) 

call x3name (' ',1) 

call y3name ('Time (msec)',11) 

call z3name ('Amplitude',9) 

call vuangl (0.,0. ,30. ) 

call yintax 

call zintax 

call xnonum 

call graf3d (0.,mfreq/4.,mfreq,0.,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,' Second subplot axes are complete' 
call reset ('hwshd') 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 
call endgr (3) 

print*,' Second subplot surface is complete' 
c THIRD SUBPLOT 
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call reset ('xnonum') 
call height (.150) 
call physor (.5,.625) 
call area2d (3.5,4.875) 
call blsur 

call volm3d (8. ,8. ,9. ) 

call x3name ('Frequency (Hz)',15) 

call y3name (' ',1) 

call z3name ('Amplitude',9) 

call vuangl (-90.,0.,30.) 

call xintax 

call zintax 

call ynonum 

call graf3d (0.,mfreq/4.,mfreq,0.,mtime/2.,mtime,zmin, 
&(zmax-zmin)/3.,zmax) 
print*,' Third subplot axes are complete' 
call reset ('hwshd') 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 
call endgr (4) 

print*,' Third subplot surface is complete' 

C FOURTH SUBPLOT 

call reset ('ynonum') 
call height (.125) 
call physor (4.75,1.5) 
call area2d (2.75,3.0) 
call blsur 

call xname ('Frequency (Hz)',15) 
call yname ('Time (msec)',11) 
call xintax 
call yintax 

call graf (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime) 
print*,' Fourth subplot axes are complete' 
call messag ('Contour at zmax/3',17,0.5,3.3) 
call reset ('hwshd') 
call conbgn 

call zrange (zmax/3. ,zmax/3. ) 
call conmak (rwdf,rdp2,rdp,1) 
call conend 

call conlin (1,'solid','nolabels',1,10) 
call contur (1,'nolabels','draw') 
call endgr (5) 

print*,' Fourth subplot surface is complete' 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [rossano.thesis]plot3dsplit64.include 

c 

c This subroutine uses disspla 11. 0 to plot a series of 

c 3 dimensional graphs of rwdf(128,64) 


subroutine plot3dsplit64 (rwdf,outname,mtime,mfreq,dp,tcode, 
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&title,mnq,ampl) 

real rotime,mfreq,zmax,zmin,rwdf(128,64),ampl 
character*23,datetime,tcode,dpnum,amplif 
character*50,title,outname,meani 
integer*4 status 

int eger ip lot_val, dp, niesh, rdp2, rdp, mnq 
rdp=64 

if (dp. ecj. 128) then 

dpnum- 128 data points' 
endif 

if (dp.eq.256) then 

dpnum- 256 data points’ 
endif 

if (dp.eq. 512) then 

dpnum- 512 data points' 
endif 

if (dp.eq.1024) then 

dpnum-1024 data points' 
endif 

if (mnq.eq. 1) then 

meani='Mean value removed' 
endif 

if (mnq.eq. 2) then 

meani='Mean value not removed' 
endif 

amplif='Time amplified' 

if (ampl.eq.0.0) amplif='Time amplified by 0. O' 
if (ampl.eq.1.0) amplif='Time amplified by 1.0' 
if (ampl.eq. 2. 0) amplif='Time amplified by 2.0' 
if (ampl.eq.3. 0) amplif='Time amplified by 3.0' 
if ( ampl. eq. 4. 0) amplif='Time amplified by 4.0' 
if (ampl.eq.5. 0) amplif='Time amplified by 5.0' 
call range64 (rwdf,rdp,zmax,zmin) 

print*,' Enter limits for plotting (do not forget decimal)' 

print*,' ZMIN ZMAX 

write(5,*) zmin.zmax 

read(5,*) zmin,zmax 

print* 

mesh=l 

rdp2=rdp*2 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val . eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val .eq. 2) then 

print*,' Please be patient' 

print*,' This will take several minutes.' 
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call LN03I 
end if 

call swissm 
call hwshd 
call chrpat (16) 

c LABELS 

call height (0.200) 

call physor (0.5,0.625) 

call area2d (7.5,9.75) 

call alnmes (0.5,0.0) 

call messag (title,100,7.5/2.,9.75) 

call messag (outname,100,7. 5/2. ,9. 5) 

call reset ('alnmes') 

call height (0.150) 

call alnmes (0.0,1.0) 

call messag (meani,50,5. 0,0. 0) 

call messag (amplif,23,5. 0,-0. 2) 

call messag (tcode,23,5. 0,-0. 4) 

status = lib$date_time (datetime) 

call messag (datetime,23,0. 0,0. 2) 

call messag (dpnum,23,0. 0,0. 0) 

call messag ('Reduced to 128 x 64',23,0.0,-0.2) 

call messag ('Smoothed 10 x 10',23,0.0,-0.4) 

call reset ('alnmes') 

call endgr (1) 

print*,' Labels are complete' 


c FIRST SUBPLOT 

call height (.150) 
call physor (.5,5.5) 
call area2d (3.5,4.875) 
call blsur 

call volm3d (8.,8.,9.) 

call x3name ('Frequency (Hz)’,15) 

call y3name ('Time (msec)’,11) 

call z3name ('Amplitude',9) 

call vuangl (-50. ,0. ,30. ) 

call xintax 

call yintax 

call zintax 

call graf3d (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,' First subplot axes are complete' 
call reset ('hwshd') 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 
call endgr (2) 

print*,' First subplot surface is complete' 


c SECOND SUBPLOT 
call height 
call physor 
call area2d 
call blsur 
call volm3d 
call x3name 


(•150) 

(4.25.5.5) 

(3.5.4.5) 

(8. ,8. ,9. ) 

(’ M) 
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call y3name ('Time (msec)',11) 
call z3name ('Amplitude',9) 
call vuangl (0. ,0. ,30. ) 
call yintax 
call zintax 
call xnonum 

call graf3d (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,' Second subplot axes are complete' 
call reset ('hwshd') 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 
call endgr (3) 

print*,' Second subplot surface is complete' 


c THIRD SUBPLOT 

call reset ('xnonum') 
call height (.150) 
call physor (.5,.625) 
call area2d (3.5,4.875) 
call blsur 

call volm3d (8. ,8. ,9. ) 

call x3name ('Frequency (Hz)’,15) 

call y3name (' ' , 1) 

call z3name ('Amplitude',9) 

call vuangl (-90. ,0. ,30. ) 

call xintax 

call zintax 

call ynonum 

call graf3d (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3. ,zmax) 
print*,' Third subplot axes are complete' 
call reset ('hwshd j 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 
call endgr (4) 

print*,' Third subplot surface is complete' 


C FOURTH SUBPLOT 

call reset ('ynonum') 
call height (.125) 
call physor (4.75,1.5) 
call area2d (2.75,3.0) 
call blsur 

call xname ('Frequency (Hz)',15) 
call yname ('Time (msec)',11) 
call xintax 
call yintax 

call graf (0. ,mfreq/4. ,mfreq,0.,mtime/2. ,mtime) 
print*,' Fourth subplot axes are complete' 
call messag ('Contour at zmax/3',17,0.5,3.3) 
call reset ('hwshd') 
call conbgn 

call zrange (zmax/3. ,zmax/3. ) 
call conraak (rwdf,rdp2,rdp,l) 
call conend 

call conlin (1,'solid','nolabels1,10) 
call contur (1,'nolabels','draw') 
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call endgr (5) 

print*,' Fourth subplot surface is complete* 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [rossano.thesis]plot3dsplitl28.include 

c 

c This subroutine uses disspla 11.0 to plot a series of 

c 3 dimensional graphs of rwdf(256,128) 


subroutine plot3dsplitl28 (rwdf,outname,mtime,mfreq,dp,tcode, 
&title,mnq,ampl) 

real mtime,mfreq,zmax,zmin,rwdf(256,128),ampl 
character*23,datetime,tcode,dpnum,amplif 
character*50,title,outname,meani 
integer*4 status 

integer iplot_val,dp,mesh,rdp2,rdp,mnq 
rdp=128 

if (dp. eq. 128) then 

dpnum- 128 data points' 
end if 

if (dp.e^.256) then 

dpnum- 256 data points' 
end if 

if (dp.eq. 512) then 

dpnum='512 data points' 
end if 

if (dp.eq.1024) then 

dpnum-1024 data points' 
end if 

if (mnq. eq. 1) then 

meani='Mean value removed' 
endif 

if (mnq. eq. 2) then 

meani='Mean value not removed' 
endif 

amplif='Time amplified' 

if (ampl.eq.0.0) amplif='Time amplified by 0. O' 
if (ampl.eq.1.0) amplif='Time amplified by 1.0' 
if (ampl.eq.2.0) amplif='Time amplified by 2.0' 
if (ampl.eq.3.0) amplif='Time amplified by 3.0' 
if (ampl.eq.4.0) amplif='Time amplified by 4.0' 
if (ampl.eq.5.0) amplif='Time amplified by 5.0' 
call rangel28 (rwdf,rdp,zmax,zmin) 

print*,' Enter limits for plotting (do not forget decimal)' 

print*,' ZMIN ZMAX r 

write(5,*) zmin,zmax 

read(5,*) zmin,zmax 

print* 

mesh=l 
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rdp2=rdp*2 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ’ (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val . eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val .eq. 2) then 

print*,' Please be patient' 
print*,’ This will take several minutes.' 
call LN03I 
endif 

call swissm 
call hwshd 
call chrpat (16) 

c LABELS 

call height (0.200) 

call physor (0.5,0.625) 

call area2d (7.5,9.75) 

call alnmes (0.5,0.0) 

call messag (title,100,7.5/2. ,9.75) 

call messag (outname,100,7.5/2.,9.5) 

call reset ('alnmes') 

call height (0.150) 

call alnmes (0.0,1.0) 

call messag (meani,50,5.0,0.0) 

call messag (amplif,23,5.0,-0.2) 

call messag (tcode,23,5.0,-0.4) 

status = lib$date_time (datetime) 

call messag (datetime,23,0.0,0.2) 

call messag (dpnum,23,0. 0,0.0) 

call messag ('Reduced to 256 x 128*,23,0.0,-0.2) 
call messag ('Smoothed 10 x 10',23,0.0,-0.4) 
call reset ('alnmes') 
call endgr (1) 

print*,' Labels are complete' 


c FIRST SUBPLOT 

call height 
call physor 
call area2d 
call blsur 
call volm3d 
call x3name 
call y3name 
call z3name 
call vuangl 
call xintax 
call yintax 
call zintax 


(.150) 

(.5,5.5) 

(3.5,4.875) 

(8. ,8. ,9.) 

('Frequency (Hz)',13) 
('Time (msec)',11) 

('Amplitude',9) 

(-50. ,0. ,30. ) 


135 







call graf3d (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3.,zmax) 
print*,' First subplot axes are complete' 
call reset ('hwshd j 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 
call endgr (2) 

print*,' First subplot surface is complete' 

c SECOND SUBPLOT 

call height (.150) 
call physor (4.25,5.5) 
call area2d (3.5,4.5) 
call blsur 

call volm3d (8. ,8. ,9. ) 

call x3name (' ,1) 

call y3name ('Time (msec)’,11) 

call z3name ('Amplitude',9) 

call vuangl (0. ,0. ,30. ) 

call yintax 

call zintax 

call xnonum 

call graf3d (0.,mfreq/4.,mfreq,0.,mtime/2.,mtime,zmin, 
&(zmax-zmin)/3.,zmax) 
print*,' Second subplot axes are complete' 
call reset ('hwshd') 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 
call endgr (3) 

print*,' Second subplot surface is complete' 

c THIRD SUBPLOT 

call reset ('xnonum') 
call height (.150) 
call physor (.5,.625) 
call area2d (3.5,4.875) 
call blsur 

call volm3d (8. ,8. ,9. ) 

call x3name ('Frequency (Hz)',15) 

call y3name (' ',1) 

call z3name ('Amplitude',9) 

call vuangl (-90.,0.,30.) 

call xintax 

call zintax 

call ynonum 

call graf3d (0. ,mfreq/4. ,mfreq,0.,mtime/2. ,mtime,zmin, 
&(zmax-zmin)/3.,zmax) 
print*,' Third subplot axes are complete' 
call reset ('hwshd') 

call surmat (rwdf,mesh,rdp2,mesh,rdp,0) 
call endgr (4) 

print*,' Third subplot surface is complete' 

C FOURTH SUBPLOT 

call reset ('ynonum') 
call height (.125) 
call physor (4.75,1.5) 
call area2d (2.75,3.0) 
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call blsur 

call xname ('Frequency (Hz)’,15) 
call yname ('Time (msec)',11) 
call xintax 
call yintax 

call graf (0.,mfreq/4. ,mfreq,0.,mtime/2. ,mtime) 
print*,' Fourth subplot axes are complete' 
call messag ('Contour at zmax/3',17,0.5,3.3) 
call reset ('hwshd') 
call conbgn 

call zrange (zmax/3. ,zmax/3. ) 
call conmak (rwdf,rdp2,rdp,1) 
call conend 

call conlin (1,'solid','nolabels',1,10) 
call contur (1,'nolabels','draw') 
call endgr (5) 

print*,' Fourth subplot surface is complete' 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [rossano. thesis]plotcon32. include 

c 

c This subroutine uses flisspla 11.0 to plot contour 

c graphs of rwdf(64,32) 


subroutine plotcon32 (rwdf,outname,mtime,mfreq,dp,tcode, 
&title,mnq,ampl) 

real mtime,mfreq,zmax,zmin,rwdf(64,32),ampl 
character*23,datetime,tcode,dpnum,amplif 
character*50,title,outname,meani 
integer*4 status 

integer iplot_val,dp,mesh,rdp2,rdp,mnq 
rdp=32 

if (dp.eq.128) then 

dpnum='l28 data points' 
endif 

if (dp.eq.256) then 

dpnum=’256 data points' 
endif 

if (dp.eq. 512) then 

dpnum='512 data points' 
endif 

if (dp.eq. 1024) then 

dpnum='l024 data points' 
endif 

if (mnq.eq. 1) then 

meani='Mean value removed' 
endif 

if (mnq.eq.2) then 

meani='.M€an value not removed' 
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endif 

amplify Time amplified' 

if (ampl.eq.0.0) amplif='Time amplified by 0.O' 
if (ampl.eq. 1. 0) amplif='Time amplified by 1.0' 
if (ampl. eq. 2. 0) amplif='Time amplified by 2.0' 
if (ampl.eq.3.0) amplif='Time amplified by 3.0' 
if (ampl.eq.4.0) amplif='Time amplified by 4.0' 
if (ampl.eq.5.0) amplif='Time amplified by 5.0' 
call range32 (rwdf,rdp,zmax,zmin) 

print*,' Enter limits for plotting (do not forget decimal)' 

print*,' ZMIN ZMAX r 

write(5,*) zmin,zmax 

read(5,*) zmin,zmax 

print* 

mesh=T 

rdp2=rdp*2 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val . eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val .eq. 2) then 

print*,' Please be patient' 
print*,' This will take several minutes.' 
call LN03I 
endif 

call swissm 
call hwshd 
call chrpat (16) 

c LABELS 

call height (0.200) 

call physor (0.5,0.625) 

call area2d (7.5,9.75) 

call alnmes (0.5,0.0) 

call messag (title,100,7.5/2. ,9.75) 

call messag (outname,100,7.5/2.,9.5) 

call messag ('Contours from zmax/6 to zmax',100,7.5/2.,9.25) 

call reset ('alnmes') 

call height (0.150) 

call alnmes (0.0,1.0) 

call messag (meani,50,5.0,0.0) 

call messag (amplif,23,5.0,-0.2) 

call messag (tcode,23,5.0,-0.4) 

status = lib$date_time (datetime) 

call messag (datetime,23,0.0,0.2) 

call messag (dpnum,23,0.0,0.0) 

call messag ('Reduced to 64 x 32',23,0.0,-0.2) 

call messag ('Smoothed 10 x 10',23,0. 0,-0.4) 

call reset ('alnmes') 
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call endgr (1) 

print*,' Labels are complete' 


c PLOT 

call height (.125) 
call physor (0.75,1.5) 
call area2d (7.0,8. G) 
call blsur 

call xname ('Frequency (Hz)',15) 

call yname ('Time (msec)',11) 

call xintax 

call yintax 

call xticks (10) 

call yticks (10) 

call graf (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime) 
print*,' Axes are complete' 
call reset ('hwshd') 
call conbgn 

call zrange (zmax/6. ,zmax) 

call conmak (rwdf,rdp2,rdp,zmax/6. ) 

call conend 

call conlin (1,'solid','nolabels',1,10) 
call conlin (2,'dot','nolabels',1,10) 
call conlin (3,'solid','nolabels’,1,10) 
call conlin (4,'dot','nolabels',1,10) 
call conlin (5,'solid','nolabels',1,10) 
call conlin (6,'dot','nolabels'^1,10) 
call contur (6,'nolabels','draw') 
call endgr (2) 

print*,' Contours are complete' 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [rossano.thesis]plotcon64. include 

c 

c This subroutine uses disspla 11.0 to plot contour 

c graphs of rwdf(128,64) 


subroutine plotcon64 (rwdf,outname,mtime,mfreq,dp,tcode, 
&title,mnq,ampl) 

real mtime,mfreq,zmax,zmin,rwdf(128,64),ampl 
character*23,datetime,tcode,dpnum,amplif 
character*50,title,outname,meani 
integer*4 status 

integer iplot_val,dp,mesh,rdp2,rdp,mnq 
rdp=64 

if (dp. eq. 128) then 

dpnum- 128 data points’ 
endif 

if (dp. eq. 256) then 
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dpnum='256 data points' 
end if 

if (dp.eq.512) then 

dpnum- 512 data points' 
endif 

if (dp. ecp 1024) then 

dpnum-1024 data points' 
endif 

if (mnq. eq. 1) then 

meani='Mean value removed' 
endif 

if (mnq.eq. 2) then 

meani='Mean value not removed' 
endif 

smplif='Time amplified' 
if (ampl.eq. 0. 0) amplif= 
if (ampl.eq. 1. 0) amplif= 
if (ampl.eq. 2. 0) amplif= 
if (ampl. eq. 3. 0) amplif= 
if (ampl. eq. 4. 0) amplif= 
if (ampl. eq. 5. 0) amplif= 


Time amplified by 0. O' 
Time amplified by 1.0' 
amplified by 2. O' 
amplified by 3. O' 
amplified by 4. O' 
amplified by 5.0' 


Time 

Time 

Time 

Time 


call range64 (rwdf,rdp,zmax,zmin) 
print*,' Enter limits for plotting (do not forget decimal)' 
™ TV ZMAX 


print*,' ZMIN 

write(5,*) zmin,zmax 
read(5,*) zmin,zmax 
print* 
mesh=l 


rdp2=rdp*2 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val .eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val .eq. 2) then 

print*,' Please be patient' 
print*,' This will take several minutes.' 
call LN03I 
endif 

call swissm 
call hwshd 
call chrpat (16) 


get 


c LABELS 

call 

call 

call 

call 

call 

call 


height 

(0. 200) 


physor 

(0. 5,0. 625) 


area2d 

(7.5,9. 75) 


alnmes 

(0. 5,0. 0) 


messag 

(title,100,7.5/2. ,9. 

75) 

messag 

(outname,100,7. 5/2. , 

9.5) 
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call messag ('Contours from zmax/6 to zmax',100,7. 5/2.,9. 25) 

call reset ('alnmes') 

call height (0.150) 

call alnmes (0. 0,1.0) 

call messag (meani,50,5.0,0.0) 

call messag (amplif,23,5. 0,-0. 2) 

call messag (tcode,23,5.0,-0.4) 

status = lib$date_time (datetime) 

call messag (datetime,23,0. 0,0. 2) 

call messag (dpnum,23,0.0,0. 0) 

call messag ('Reduced to 128 x 64',23,0.0,-0.2) 

call messag ('Smoothed 10 x 10',23,0.0,-0.4) 

call reset ('alnmes') 

call endgr (1) 

print*,' Labels are complete' 


c PLOT 

call height (.125) 
call physor (0.75,1.5) 
call area2d (7.0,8.0) 
call blsur 

call xname ('Frequency (Hz)',15) 

call yname ('Time (msec)',11) 

call xintax 

call yintax 

call xticks (10) 

call yticks (10) 

call graf (0. ,mfreq/4. ,mfreq,0. ,mtime/2. ,mtime) 
print*,' Axes are complete' 
call reset ('hwshd') 
call conbgn 

call zrange tzmax/6.,zmax) 

call conmak (rwdf,rdp2,rdp,zmax/6. ) 

call conend 

call conlin (1,'solid','nolabels',1,10) 
call conlin (2,'dot','nolabels',1,10) 
call conlin (3,'solid','nolabels',1,10) 
call conlin (4,'dot','nolabels',1,10) 
call conlin (5,'solid','nolabels',1,10) 
call conlin (6,'dot','nolabels',1,10) 
call contur (6,'nolabels','draw') 
call endgr (2) 

print*,' Contours are complete' 
call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [rossano. thesis]plotconl28. include 

c 

c This subroutine uses disspla 11.0 to plot contour 

c graphs of rwdf(256,128) 


subroutine plotconl28 (rwdf,outname,mtime,mfreq,dp,tcode, 
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&title,mnq,amp1) 

real mtime,mfreq,zmax,zmin,rwdf(256,128),ampl 
character*23,datetime,tcode,dpnum,amplif 
character*50,title,outname,meani 
integer*4 status 

integer iplot_val,dp,mesh,rdp2,rdp,mnq 
rdp=128 

if (dp. ecj. 128) then 

dpnum- 128 data points' 
end if 

if (dp. ecj. 256) then 

dpnum= 256 data points' 
endif 

if (dp. eq. 512) then 

dpnum- 512 data points' 
endif 

if (dp.eq.1024) then 

dpnum- 1024 data points' 
endif 

if (mnq.eq.1) then 

meani='Mean value removed' 
endif 

if (mnq. eq. 2) then 

meani='Mean value not removed' 
endif 

amplif='Time amplified’ 

if (ampl.eq.0.0) amplif='Time amplified by 0. O' 
if (ampl.eq.1.0) amplif='Time amplified by 1.0' 
if (ampl.eq.2.0) amplif='Time amplified by 2.0' 
if (ampl.eq.3. 0) amplif='Time amplified by 3.0' 
if (ampl.eq.4.0) amplif='Time amplified by 4.0' 
if (ampl.eq.5.0) amplif='Time amplified by 5.0’ 
call rangel28 (rwdf,rdp,zmax,zmin) 

print*,' Enter limits for plotting (do not forget decimal)' 

print*,’ ZMIN ZMAX' 

write(5,*) zmin,zmax 

read(5,*) zmin,zmax 

print* 

mesh=l 

rdp2=rdp*2 

call reset ('all') 

write(5,*) ' Do you want to view the plot on the screen or get 
& a hardcopy?' 

write(5,*) ' (1 for view, 2 for hardcopy)' 

print* 

read(5,*) iplot_val 
print* 

if (iplot_val . eq. 1) then 

print*,' When finished viewing hit return key' 
call pgpx 
endif 

if (iplot_val .eq. 2) then 

print*,' Please be patient' 

print*,' This will take several minutes.' 
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call LN03I 
endii 

call swissm 
call hwshd 
call chrpat (16) 


c LABELS 

call height (0.200) 

call physor (0.5,0.625) 

call area2d (7.5,9.75) 

call alnmes (0.5,0.0) 

call messag (title,100,7. 5/2. ,9.75) 

call messag (outname,100,7. 5/2. ,9.5) 

call messag ('Contours from zmax/6 to zmax', 100,7. 5/2. 

call reset ('alnmes') 

call height (0.150) 

call alnmes (0. 0,1.0) 

call messag (meani,50,5. 0,0. 0) 

call messag (amplif,23,5. 0,-0. 2) 

call messag (tcode,23,5. 0,-0. 4) 

status = lib$date_time (datetime) 

call messag (datetime,23,0. 0,0. 2) 

call messag (dpnum,23,0.0,0. 0) 

call messag ('Reduced to 256 x 128',23,0.0,-0.2) 
call messag ('Smoothed 10 x 10*,23,0.0,-0.4) 
call reset ('alnmes') 
call endgr (1) 

print*,' Labels are complete' 


c PLOT 

call height (.125) 
call physor (0.75,1.5) 
ca)l area2d (7.0,8. 0) 
call blsur 

call xname ('Frequency (Hz)',15) 

call yname ('Time (msec)',11) 

call xintax 

call yintax 

call xticks (10) 

call yticks (10) 

call graf (0. ,mfreq/4. ,mfreq,0. ,ratime/2. ,mtime) 
print*,' Axes are complete' 
call reset ('hwshd') 
call conbgn 

call zrange (zmax/6.,zmax) 

call conmak (rwdf,rdp2,rdp,zmax/6. ) 

call conend 

call conlin (1,'solid','nolabels',1,10) 
call conlin (2,'dot','nolabels',1,10) 
call conlin (3,'solid','nolabels',1,10) 
call conlin (4,'dot',’nolabels',1,10) 
call conlin (5,'solid','nolabels',1,10) 
call conlin (6,'dot','nolabels',1,10) 
call contur (6,'nolabels','draw') 
call endgr (2) 

print*,' Contours are complete' 


.25) 
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call endpl (0) 
print*,' Plotting complete' 
print* 
return 
end 
c 

c [rossano.thesis]pseudo, include 

c 

c This subroutine smooths the WDF along both the time and 

c frequency axes 

c 


subroutine pseudo (dimt,dimf,rwdf,dp,dt,nn,mm) 

integer dp,i,j,dimt,dimf,dp2,rdp,rdp2,nn,mm,nf,mt,nf2,mt2, 
&ii,jj,L,LL,k,kk 

real rwdf(256,128),hg(-25: 25,-25: 25),pi,dt,df, 

&wdf(1100,550),coef.val 
common /wdfc/ wdf 

pi=4.*atan(1. ) 
df=l. /(4. *dp*dt) 
dp2=dp*2 
rdp=dp/mm 
rdp2=dp2/nn 

print*,' Smoothing 10 x 10' 

nf=10 

mt=10 

nf2=nf*2 

mt2=mt*2 

val=l. /(C2. *pi*nf*df)*(2.*pi*mt*dt)) 
do 20 j=-mt2,mt2 
do 10 i=-nf2,nf2 

coef=-((j*j)/(2. *mt*mt))-((i*i)/(2. *nf*nf)) 
hg(i,j)=val*exp(coef) 


10 

continue 

20 

continue 


do 100 j=l,rdp 


do 100 i=l,rdp2 

100 

rwdf(i,j)=0. 0 


do 500 i=l,dp2,nn 


ii=l+(i-l)/nn 


do 450 j=l,dp,mm 
jj=l+(j-l)/mm 
do 400 L=i-nf2,i+nf2 
LL=L 

if (L. It. 1) LL=l+dp2 
if (L. gt.dp2) LIr=l-dp2 
do 350 k=j-rat2,j+mt2 
kk=k 

if (k.It.1) kk=k+dp 
if (k.gt.dp) kk=k-dp 

rwdf(ii,jj)=rwdf(ii,jj)+wdf(LL,kk)*hg(L-i,k-j) 
350 continue 

400 continue 

450 continue 
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500 continue 

c This loop removes a plotting artifact 
do 600 j=l,rdp 

rwdf(rdp2-3,j)=0. 0 
rwdf(rdp2-2,j)=0. 0 
rwdf(rdp2-l,j)=0. 0 
rwdf(rdp2,j)=0. 0 
600 continue 
return 
end 
c 

c [rossano.thesis]range32. include 

c 

c This subroutine finds the maximum and minimum amplitudes 

c of array rwdf(64,32) 


subroutine range32 (rwdf,dp,zmax,zmin) 

real zmax,zmin,rwdf(64,32) 
integer i,j,dp,dp2 

dp2=dp*2 
zmax=0. 0 
zmin=0. 0 
do 100 j=l,dp 
do 200 i=l,dp2 

if (rwdf(i,j).gt.zmax) zmax=rwdf(i,j) 
if (rwdf(i, j). It. zmin) znin=rwdf(i,j) 

200 continue 
100 continue 
return 
end 
c 

c [ rossano. thesis]range64. include 

c 

c This subroutine fin^-. the maximum and minimum amplitudes 

c of array rwdf(128,6 


subroutine range64 (rwdf,dp,zmax,zmin) 

real zmax,zmin,rwdf(128,64) 
integer i,j,dp,dp2 

dp2=dp*2 
zmax=0. 0 
zmin=0.0 
do 100 j=l,dp 
do 200 i=l,dp2 

if (rwdf(i,j). gt. zmax) zmax=rwdf(i,j) 
if (rwdf(i,j). It. zmin) zmin=rwdf(i,j) 
200 continue 
100 continue 
return 
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end 


c [ rossano. thesis]rangel28. include 

c 

c This subroutine finds the maximum and minimum amplitudes 

c of array rwdf(256,128) 


subroutine rangel28 (rwdf,dp,zmax,zmin) 

real zmax,zmin,rwdf(256,128) 
integer i,j,dp,dp2 

dp2=dp*2 
zmax=0. 0 
zmin=0. 0 
do 100 j=l,dp 
do 200 i=l,dp2 

if (rwdf(i, j). gt. zrnax) zmax=rwdf(i,j) 
if (rwdf(i,j). It. zmin) zrain=rwdf(i,j) 

200 continue 

100 continue 
return 
end 
c 

c [rossano.thesis]reduce.include 

c 

c This subroutine reduces the data in wdf to rwdf 


subroutine reduce (dimt,dimf,dt,dp,nn,mm,rwdf,rval) 

integer dimt,dimf,dp,dp2,i,j,ii,jj,nn,mm,rval 
real wdf(1100,550),rwdf(256,128) 
common /wdfc/ wdf 

dp2=dp*2 
df=l. /(4. *dp*dt) 


print*,' 

What do 

you 

want to reduce to?' 

print*,' 

1 . 

64 

x 32' 

print*,' 

2. 

128 

x 64' 

print*,' 

3. 

256 

x 128' 


read(5,*) rval 

if (rval.eq. 1) then 
nn=dp2/64 
mm=dp/32 
endif 

if (rval.eq. 2) then 
nn=dp2/128 
mm=dp/64 
endif 

if (rval.eq. 3) then 
nn=dp2/256 
mm=dp/128 
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endif 
ii=0 
j j=0 

do 100 j=1,dp,mm 
ii=0 

do 200 i=l,dp2,nn 
ii=ii+l 

rwdf(ii,jj)=wdf(i,j) 

200 continue 
100 continue 
return 
end 
c 

c [ rossano.thesis]size, include 

c 

c This code is included in wigfunl.for and wigfunlb. for and their 

c subroutines. It provides an easy way to change the size of the 

c time plotting arrays, 

c 

c real x(128),y(128),t(128) 

c real x(256),y(256),t(256) 

real x(512),y(512),t(512) 
c real x(1024),y(1024),t(1024) 

integer sizedp 

sizedp=512 

c 

c [rossano. thesis]timesig. include 

c 

c This subroutine modifies and plots the time signal 

c 


subroutine timesig (dimt,ain,dp,dt,outname,tcode, 
&title,s,sr,si,mnq,amp1) 

integer dp,mnq,wndwcode,anq,dimt,atvq 
real ain(dimt),sr(dimt),si(dirat),meanv,dt,zmax,ampl 
complex s(dimt) 
character*23 tcode 

character*50 title,titlehold,outname 
print* 

print*,' Do you want to plot the raw time signal?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
if (atvq. eq. 1) then 

title='Raw Time Signal' 
call plot2d (ain,dt,dp,title,outname) 
endif 
print* 

call mean (dimt,ain,dp,meanv) 
print*,' The mean value is' 
write(5,906) meanv 
print* 
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print*,' Do you want to remove the mean from the time signal?' 
print*,’ (1 for yes, 2 for no)' 

read(*,*) mnq 

if (mnq. eq.1) call meanr (dimt,ain,dp,meanv) 

call maxamp (ain,dp,zmax) 

print* 

print*,' The max amplitude is' 

print*,zmax 

print* 

print*,' You want this to be approx 1 in order to avoid 
&plotting artifacts. ' 
ampl=l. 0/zmax 
print* 

print*,' Recommend an amplification of' 

print*,ampl 

print* 

print*,' Do you want to amplify the signal?' 
print*,' (1 for yes, 2 for no)’ 

read(5,*) atvq 

if (atvq.eq.1) call amplify (dimt,ain,dp,ampl) 

if (atvq. ne. 1) amp 1=0. 0 

print* 

print*,' Do you want to window the time signal?' 
print*,' (1 for yes, 2 for no) r 

read(5,*) atvq 
if (atvq. eq.1) then 

call window (dimt,ain,dp,dt,wndwcode) 
if (wndwcode.eq.2) tcode='Hanning window time' 
if (wndwcode. eq. 3) tcode='Hamming window’ time' 
endif 
print* 

print*,' Do you want to plot the modified time signal?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
if (atvq.eq.1) then 

title='Modified Time Signal' 
call plot2d (ain.dt,dp,title,outname) 
endif 

title='Wigner Distribution' 
print* 

print*,' Do you want to make the time signal analytic?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) anq 
if (anq.eq.1) then 

title='Wigner-Ville Distribution' 
call analytic (dimt,ain,s,dp) 
do 100 j=l,dp 
sr(j)=ain(j) 
si(j)=aimag(s(j)) 

100 continue 
print* 

print*,' Do you want to plot the analytic time signal?' 
print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
if (atvq.eq.1) then 
print* 
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if (dp. ge.512) then 

print*,' The STD00001.dat file may be too large to print.' 
print*,' For a hardcopy run [rossano. thesis]wigfunlb. for' 
print* 

print*,' Do you want to continue with analytic plotting 

& here?' 

print*,' (1 for yes, 2 for no)' 

read(5,*) atvq 
if (atvq.eq.2) go to 300 
endif 

titlehold=title 
title='Analytic Time Signal' 
call plot2d2 (sr,si,dt,dp,title,outname,mnq) 
title=titlehold 
endif 
endif 

300 if (anq.eq.2) then 
do 200 j=l,dp 
sr(j)=ain(j) 
si( j)=0. 0 

s(j)=cmplx(sr(j),si(j)) 


200 

continue 


endif 


return 

906 

format(f16.8) 


end 


c [rossano. thesis]wigh. include 

c 

c This subroutine calculates the WDF 


subroutine wigh (dimt,dimf,s,c,df,dt,dp2,dp) 

integer dimt,dimf,i,j,dp2,dp 
complex s(dimt),c(dimf) 
real wdf(1100,550),df,dt 
common /wdfc/ wdf 

do 100 j=l,dp 

call corr (dimt,dimf,s,j,dt,c,dp) 
call fft (dimf,c,dp2,0) 
do 200 i=l,dp2 

wdf(i,j)=real(c(i)) 

200 continue 

100 continue 
return 
end 
c 

c [rossano.thesis]window.include 

c 

c This subroutine calls the available windowing functions 


subroutine window (dimt,ain,dp,dt,wndwcode) 
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integer atvq,dp,wndwcode,dimt 
real ain(dimt),dt 

print* 

print*,' Which window would you like to apply?' 

print*,' 1. None' 

print*,' 2. Hanning' 

print*,' 3. Hamming' 

read(5,*) atvq 

if (atvq.eq. 1) go to 100 

if (atvq.eq.2) then 

call hanning (dimt,ain,dp,dt) 
wndwcode=2 
end if 

if (atvq.eq. 3) then 

call hamming (dimt,ain,dp,dt) 
wndwcode=3 
end if 

100 return 
end 

D. DATA CONVERSION PROGRAM 

c 

c [rossano.thesis]convert, for 

c 

c This program converts the experimental input data file obtained 

c from HP Vista into the format used for WIGFUN1. FOR 


real tim(2),amp(2) 
integer modq,i,n 
character*20 inname,outname 
character*l c(6) 

print* 

print*,' This program converts HP Vista data input format into' 
print*,' the data format used in 
& [rossano.thesis]VIGFUN1. FOR' 
print* 

print*,' Enter HP Vista input filename' 

read(5,901) inname 

print* 

print*,' Has the last line been modified so that 9.999999' 
print*,' is in each column to indicate an end of file?' 
print*,' An example follows:' 
print*,' ( 4.960938e-01, 7.668274e-01) 

& ( 4.980469e-01, 1.742336e-01)' 

print*,' ( 5.000000e-01, -5.128824e-01) 

& ( 9.999999 , 9.999999 )' 

print*,’ ( 9.999999 , 9.999999 ) 

& ( 9.999999 , 9.999999 )' 

print*,' (1 for yes, 2 for no)' 

read(5,902)modq 
if (modq.eq. 2) then 

print*,' You need to edit the HP Vista file’ 
print*,’ Good bye' 
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go to 900 
endif 
print* 

print*,’ Enter output filename (i.e. sin62.5e)' 
read(5,901) outname 

open (unit=4,file=inname,status='old') 
rewind 4 

open (unit=7,file=outname,status='new') 
print* 

print*,' Delta t = ' 

print* 

n=0 

do 100 i=l,10000 

read(4,903) c(1),tim(1),c(2),amp(1),c(3),c(4),tim(2), 

& c(5),amp(2),c(6) 

if (tim(1). eq.9.999999) then 

if (amp(1).eq.9.999999) go to 200 
endif 
n=n+l 

write(7,904) tim(1),amp(1) 
if (tim(2).eq. 9. 999999) then 

if (amp(2).eq. 9. 999999) go to 200 
endif 
n=n+l 

write(7,904) tim(2),amp(2) 
if (n.eq.2) go to 300 
if (n.eq.102) go to 300 
if (n. eq. 202) go to 300 
if (n.eq.302) go to 300 
if (n.eq.402) go to 300 
if (n.eq.502) go to 300 
if (n.eq.602) go to 300 
if (n.eq.702) go to 300 
if (n.eq.802) go to 300 
if (n.eq.902) go to 300 

100 continue 
print* 

print*,' There are more points remaining in HP Vista data file' 
print*,' increase the loop size' 

200 write(7,904) 9999. ,9999. 
print* 

print*,' Total number of data points is' 
write(5,905) n 
go to 900 

300 dt=tim(2)-tim(1) 
print*,dt 
go to 100 

900 close (unit=4) 
close (unit=7) 

901 format (a20) 

902 format (il) 

903 format (3x,al,el5.6,al,el5.6,al,2x,al,el5.6,al,el5.6,al) 

904 format (2x,el6.8,5x,el6.8) 

905 format (lx,il0) 
end 
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APPENDIX B. DATA TRANSFER FROM SOURCE TO VAX 

COMPUTER 


A. DATA SOURCE TO HP3565 COMPUTER 

Numerous sources of data were used. Analog sources included function generators 
and FM tape recordings of accelerometer, velocity, and pressure transducers. Another 
source of data resulted from direct measurements of vibrations of machinery in the lab¬ 
oratory. All signals were fed into an input module of an HP3565 computer based on the 
HP9000 350 CPU which included a dynamic signal analyzer. Hewlett-Packard s 
HP-Vista software was used to control the hardware. The input module was calibrated 
and ranged for the applicable voltages of the input signal. The digital filters were set to 
allow the frequencies of interest to pass. The primary purpose of the HP3565 computer 
system was to digitize and filter the data. 

B. HP3565 COMPUTER TO PC 

Once the data was digitized, a copy of the raw data was printed out for use in veri¬ 
fying the data upon arrival in the VAX computer. The digitized data was transferred 
from the HP3565 to a personal computer (PC) as a printed ASCII file output. The ac¬ 
tual data transfer took place over an RS232 cable connecting the laser printer port on 
the HP3565 to the PC communications port. 

C. PC TO VAX 

At the PC the data file was transferred to a 5.25" floppy disc. This data was down¬ 
loaded to the VAX computer in the Mechanical Engineering CAD CAE Lab via the 
terminal available there. 

D. EDITING DATA FILE UPON ARRIVAL IN THE VAX 

The data file which arrived in the VAX included a banner page and file header in¬ 
formation which had to be manually removed by editing. In addition, there were form¬ 
feed characters throughout the file which had to be removed. The data in the VAX was 
compared with the data printouts made by the HP3565. On the average there were 3 
or 4 corrections to a 1024 data point file due to the high data transfer rate across the 
RS232 cable (9600 Baud). A slower data transfer rate could have been used to improve 
the accuracy, but the transfer time would have taken much longer. A final line of 
characters was added to the data file to serve as an end of file indication. A computer 
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program, CONVERT.FOR (the last program in Appendix A), was used to convert the 
data from the edited IIP3565 format into the format used in subroutine 
DA TAIN.INCLUDE. The data format conversion program also provided a final check 
to ensure that all data points had been included by counting the number of data points. 
A block size setting of 1024 real or S192 real on the HP3565 resulted in 1024 or 8192 
data points, respectively. 
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